AliPhysics  48852ec (48852ec)
 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), fMultESDTPCdif(15000), 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), fMultESDTPCdif(15000), 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  if(fPileUp) {
431  AliMultSelection *multSelection = 0x0;
432  multSelection = static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
433  if(multSelection) {
434  fCentrality = multSelection->GetMultiplicityPercentile(fCentralityEstimator.Data());
435  if(fCentrality > fCentralityMin && fCentrality < fCentralityMax && multSelection->GetMultiplicityPercentile("CL1") < 90) {
438  }
440  return kTRUE;
441  } else {
443  return kFALSE;
444  }
445  }
446  }
447  else fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
448  Double_t cenB(-999);
449  // if a second centrality estimator is requited, set it
450  (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
451  if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
452  if(fQA) fCentralityNoPass->Fill(fCentrality) ;
453  return kFALSE;
454  }
455  const Int_t nGoodTracks = event->GetNumberOfTracks();
456  if(fCentralityCut2010) { // cut on outliers
457  Float_t multTPC(0.); // tpc mult estimate
458  Float_t multGlob(0.); // global multiplicity
459  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
460  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
461  if(!trackAOD) AliFatal("Not a standard AOD");
462  if (!trackAOD) continue;
463  if (!(trackAOD->TestFilterBit(1))) continue;
464  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;
465  multTPC++;
466  }
467  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
468  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
469  if(!trackAOD) AliFatal("Not a standard AOD");
470  if (!trackAOD) continue;
471  if (!(trackAOD->TestFilterBit(16))) continue;
472  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;
473  Double_t b[2] = {-99., -99.};
474  Double_t bCov[3] = {-99., -99., -99.};
475  AliAODTrack copy(*trackAOD);
476  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
477  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
478  multGlob++;
479  } //track loop
480  // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
481  if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
482  if(fQA) {
483  fMultCorAfterCuts->Fill(multGlob, multTPC);
484  fMultvsCentr->Fill(fCentrality, multTPC);
485  }
486  }
487 
488  if(fCentralityCut2011) { // cut on outliers
489  Float_t multTPC(0.); // tpc mult estimate
490  Float_t multGlob(0.); // global multiplicity
491  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
492  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
493  if(!trackAOD) AliFatal("Not a standard AOD");
494  if (!trackAOD) continue;
495  if (!(trackAOD->TestFilterBit(1))) continue;
496  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;
497  multTPC++;
498  }
499  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
500  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
501  if(!trackAOD) AliFatal("Not a standard AOD");
502  if (!trackAOD) continue;
503  if (!(trackAOD->TestFilterBit(16))) continue;
504  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;
505  Double_t b[2] = {-99., -99.};
506  Double_t bCov[3] = {-99., -99., -99.};
507  AliAODTrack copy(*trackAOD);
508  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
509  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
510  multGlob++;
511  } //track loop
512  //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
513  if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
514  if(fQA) {
515  fMultCorAfterCuts->Fill(multGlob, multTPC);
516  fMultvsCentr->Fill(fCentrality, multTPC);
517  }
518  }
519 
521  return kTRUE;
522 }
523 //_____________________________________________________________________________
525 {
526  // Initialize the Bayesian PID object for AOD
527  if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
529 }
530 //_____________________________________________________________________________
531 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
532 {
533  // Check if the particle passes the TPC TOF bayesian cut.
534  if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
536  if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
537  Float_t *probabilities = fBayesianResponse->GetProb();
538  if (probabilities[3] > fPIDConfig[6]) {
539  if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
540  return kTRUE;
541  }
542  return kFALSE;
543 }
544 //_____________________________________________________________________________
546 {
547  // check if track passes dca cut according to dca routine
548  // setup the routine as follows:
549  // fDCAConfig[0] < -1 no pt dependence
550  // fDCAConfig[0] = 0 do nothing
551  // fDCAConfig[0] > 1 pt dependent dca cut
552 // if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
553 
554  AliAODTrack copy(*track);
555  if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
556  if(fQA) {
557  Double_t b[2] = { -99., -99.};
558  Double_t bCov[3] = { -99., -99., -99.};
559  if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
560  fDCAXYQA->Fill(b[0]);
561  fDCAZQA->Fill(b[1]);
562  fDCAPrim->Fill(track->Pt(), b[0]);
563  }
564  }
565  return kTRUE;
566  }
567  Double_t b[2] = { -99., -99.};
568  Double_t bCov[3] = { -99., -99., -99.};
569  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
570  if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
571  if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
572  if(fDCAConfig[0] > .9) {
573  if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
574  Double_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
575  if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
576  if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
577  }
578  if(fQA) {
579  fDCAXYQA->Fill(b[0]);
580  fDCAZQA->Fill(b[1]);
581  fDCAPrim->Fill(track->Pt(), b[0]);
582  }
583  return kTRUE;
584 }
585 //_____________________________________________________________________________
586 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
587 {
588  // Kaon identification routine, based on multiple detectors and approaches
589 // if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
590 
591  /* if(fUsePidResponse) {
592  Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
593  fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
594  if(prob[3] > fPIDConfig[6]) return kTRUE;
595  }
596  if(!PassesDCACut(track)) return kFALSE;
597  if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
598  if(track->Pt() < fPIDConfig[1]) {
599  if(fDebug>1) cout << " ITS received track with p_t " << track->Pt() << endl;
600  // if tpc control is disabled, pure its pid
601  if(fPIDConfig[2] < 0.) {
602  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
603  return kFALSE;
604  }
605  // else, switch to ITS pid with TPC rejection of protons and pions
606  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
607  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
608  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) {
609  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
610  return kTRUE;
611  }
612  return kFALSE;
613  }
614  if((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4])) {
615  if(fDebug>1) cout << " TPC received track with p_t " << track->Pt() << endl;
616  // if its control is disabled, pure tpc pid
617  if(fPIDConfig[5] < 0.) {
618  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
619  return kFALSE;
620  }
621  // else, switch to TPC pid with ITS rejection of protons and pions
622  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
623  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
624  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) {
625  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
626  return kTRUE;
627  }
628  return kFALSE;
629  }
630  if(fDebug>1) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
631  // switch to bayesian PID
632  if (PassesTPCbayesianCut(track)) {
633  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal());fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
634  return kTRUE;
635  }*/
636  if(fQA) {
637  Float_t length = track->GetIntegratedLength();
638  Float_t time = track->GetTOFsignal();
639  Double_t beta = -.05;
640 
641  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
642 
643  fNOPID->Fill(track->P(), track->GetTPCsignal(), fCentralityWeight);
644  fNOPIDTOF->Fill(track->P(), beta, fCentralityWeight);
645  }
646 
647 
648 
649  //kaon stuff
650  Double_t nSigKTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
651  Double_t nSigKTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
652 
653  Double_t nSigmaK = 999;
654 
655  if(track->Pt() > .4) {
656  nSigmaK = TMath::Sqrt(nSigKTPC*nSigKTPC+nSigKTOF*nSigKTOF);
657  } else {
658  if(track->GetTPCsignal() > 110) nSigmaK = TMath::Abs(nSigKTPC);
659  }
660  // pion
661  Double_t nSigPiTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
662  Double_t nSigPiTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
663 
664  Double_t nSigmaPi = 999;
665 
666  if(track->Pt() > .4) {
667  nSigmaPi = TMath::Sqrt(nSigPiTPC*nSigPiTPC+nSigPiTOF*nSigPiTOF);
668  } else {
669  if(track->GetTPCsignal() < 60) nSigmaPi = TMath::Abs(nSigPiTPC);
670  }
671 
672  //proton
673  Double_t nSigPTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
674  Double_t nSigPTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
675 
676  Double_t nSigmaP = 999;
677 
678  if(track->Pt() > .4) {
679  nSigmaP = TMath::Sqrt(nSigPTPC*nSigPTPC+nSigPTOF*nSigPTOF);
680  } else {
681  if(track->GetTPCsignal() > 110) nSigmaP = TMath::Abs(nSigPTPC);
682  }
683 
684  Short_t minSigma = FindMinNSigma(nSigmaPi, nSigmaK, nSigmaP);
685 
686 
687 
688  if(minSigma == 0) return kFALSE;
689  if((nSigmaK == nSigmaPi) && ( nSigmaK == nSigmaP)) return kFALSE;
690 
691  if(minSigma == 2 &&!GetDoubleCountingK(nSigmaK, minSigma)) {
692 
693  if(fQA) {
694  Float_t length = track->GetIntegratedLength();
695  Float_t time = track->GetTOFsignal();
696  Double_t beta = -.05;
697 
698  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
699 
700  // additional selection to bruteforce throw away bad betas
701  if(beta < 0.4) return kFALSE;
702 
703  fPIDk->Fill(track->P(), track->GetTPCsignal(), fCentralityWeight);
704  if(track->Pt() > .4) fPIDTOF->Fill(track->P(), beta, fCentralityWeight);
705  }
706 
707 
708  return kTRUE;
709  }
710 
711  return kFALSE;
712 }
713 //_____________________________________________________________________________
714 template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
715 {
716  // return p_t of track pair
717  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
718  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
719  TVector3 c = a + b;
720  return c.Pt();
721 }
722 //_____________________________________________________________________________
723 template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
724 {
725  // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
726  Double_t pt = PhiPt(track1, track2);
727  if (tracktype == 0) {
728  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
729  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
730  fInvMNP[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
731  if(fQA) fPtSpectra[ptbin]->Fill(pt, fCentralityWeight);
732  }
733  }
734  }
735  if (tracktype == 1) {
736  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
737  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
738  fInvMPP[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
739  }
740  }
741  }
742  if (tracktype == 2) {
743  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
744  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
745  fInvMNN[ptbin]->Fill(InvariantMass(track1, track2), fCentralityWeight);
746  }
747  }
748  }
749 }
750 //_____________________________________________________________________________
751 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
752 {
753  // check if track meets quality cuts
754  if(!track) return kFALSE;
755  return fPOICuts->IsSelected(track);
756 }
757 //_____________________________________________________________________________
758 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
759 {
760  // Set null cuts
761 // if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
762  fCutsRP->SetEvent(event, MCEvent());
764  fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
765  fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
766  fNullCuts->SetEvent(event, MCEvent());
767 }
768 //_____________________________________________________________________________
770 {
771  // Prepare flow events
772 // if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
776  fFlowEvent->DefineDeadZone(0, 0, 0, 0);
777 }
778 //_____________________________________________________________________________
780 {
781  // UserExec: called for each event. Commented where necessary
782 // if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
783  TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
784  if(fEventMixing) {
785  MixingCandidates = new TObjArray();
786  MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
787  }
788  if (!fPIDResponse) {
789  if(fDebug > 0 ) cout << " Could not get PID response " << endl;
790  return;
791  }
792  fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
793  if (fAOD) {
794  if (!EventCut(fAOD)) return; // check for event cuts
795 
796 
797  // add extra pileup cuts for high intensity runs
798  // courtesy of alex dobrin
799  if (fPileUp){
800 
801 
802 
803  // 12 11 centrality cut for pileup
804  //
805  //Centrality
806  Float_t v0Centr = -100.;
807  Float_t cl1Centr = -100.;
808  Float_t cl0Centr = -100.;
809 
810  AliMultSelection* MultSelection = 0x0;
811  MultSelection = (AliMultSelection*) fAOD->FindListObject("MultSelection");
812  if( !MultSelection) {
813  AliWarning("AliMultSelection object not found!");
814  return;
815  } else {
816  v0Centr = MultSelection->GetMultiplicityPercentile("V0M");
817  cl1Centr = MultSelection->GetMultiplicityPercentile("CL1");
818  cl0Centr = MultSelection->GetMultiplicityPercentile("CL0");
819  }
820 
821  if (v0Centr >= 90. || v0Centr < 0)
822  return;
823 
824 
825  const Int_t nTracks = fAOD->GetNumberOfTracks();
826  Int_t multEsd = ((AliAODHeader*)fAOD->GetHeader())->GetNumberOfESDTracks();
827 
828  Int_t multTrk = 0;
829  Int_t multTrkBefC = 0;
830  Int_t multTrkTOFBefC = 0;
831  Int_t multTPC = 0;
832 
833  for (Int_t it = 0; it < nTracks; it++) {
834 
835  AliAODTrack* aodTrk = (AliAODTrack*)fAOD->GetTrack(it);
836 
837  if (!aodTrk){
838  delete aodTrk;
839  continue;
840  }
841 
842  if (aodTrk->TestFilterBit(32)){
843 
844  multTrkBefC++;
845 
846  if ( TMath::Abs(aodTrk->GetTOFsignalDz()) <= 10 && aodTrk->GetTOFsignal() >= 12000 && aodTrk->GetTOFsignal() <= 25000)
847  multTrkTOFBefC++;
848 
849  if ((TMath::Abs(aodTrk->Eta()) < TMath::Abs(fCandidateMinEta)) && (aodTrk->GetTPCNcls() >= 70) && (aodTrk->Pt() >= .2) && (aodTrk->Pt() < 20))
850  multTrk++;
851 
852  }
853 
854  if (aodTrk->TestFilterBit(32))
855  multTPC++;
856 
857  }
858 
859 
860  Float_t multTPCn = multTPC;
861  Float_t multEsdn = multEsd;
862  Float_t multESDTPCDif = multEsdn - multTPCn*3.38;
863 
864 
865 
866  if (cl0Centr < fLowCut->Eval(v0Centr))
867  return;
868 
869  if (cl0Centr > fHighCut->Eval(v0Centr))
870  return;
871 
872 
873  if (multESDTPCDif > fMultESDTPCdif)
874  return;
875 
876 
877  if (multTrkTOFBefC < fMultTOFLowCut->Eval(Float_t(multTrkBefC)))
878  return;
879 
880  if (multTrkTOFBefC > fMultTOFHighCut->Eval(Float_t(multTrkBefC)))
881  return;
882 
883 
884  if (plpMV(fAOD))
885  return;
886 
887  Short_t isPileup = fAOD->IsPileupFromSPD(3);
888  if (isPileup != 0)
889  return;
890 
891  if (((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0)
892  return;
893 
894 
895  /*
896  Int_t bc2 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
897  if (bc2 != 0)
898  return;
899 
900  Int_t bc1 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
901  if (bc1 != 0)
902  return;
903  */
904 
905  // add vertexer selection
906  AliAODVertex* vtTrc = fAOD->GetPrimaryVertex();
907  AliAODVertex* vtSPD = fAOD->GetPrimaryVertexSPD();
908  if (vtTrc->GetNContributors()<2 || vtSPD->GetNContributors()<1) return; // one of vertices is missing
909  double covTrc[6],covSPD[6];
910  vtTrc->GetCovarianceMatrix(covTrc);
911  vtSPD->GetCovarianceMatrix(covSPD);
912  double dz = vtTrc->GetZ()-vtSPD->GetZ();
913  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
914  double errTrc = TMath::Sqrt(covTrc[5]);
915  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
916  if (TMath::Abs(dz)>0.2 || TMath::Abs(nsigTot)>10 || TMath::Abs(nsigTrc)>20) return; // bad vertexing
917 
918  if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fVertexRange) return;
919  fVertex = fAOD->GetPrimaryVertex()->GetZ();
920  fVertexZ->Fill(fVertex);
921 
922  fMultvsCentr->Fill(multTPCn, v0Centr);
923 
924  fMultCorAfterCuts->Fill(v0Centr,cl0Centr);
925 
926 
927  //new function for 2015 to remove incomplete events
928  if (fAOD->IsIncompleteDAQ()) return;
929  } else {
930  if(TMath::Abs(fAOD->GetPrimaryVertex()->GetZ()) > fVertexRange) return;
931  fVertex = fAOD->GetPrimaryVertex()->GetZ();
932  fVertexZ->Fill(fVertex);
933  }
934 
935 
936  // InitializeBayesianPID(fAOD); // init event objects
937 
938 
939  double qx(0); // for tpc ep
940  double qy(0);
941  SetNullCuts(fAOD);
942  PrepareFlowEvent(fAOD->GetNumberOfTracks());
943  fCandidates->SetLast(-1);
944  if(fIsMC) IsMC(); // launch mc stuff
945  if(fQA) fEventStats->Fill(0);
946  Int_t unTracks = fAOD->GetNumberOfTracks();
947  AliAODTrack* un[unTracks];
948  AliAODTrack* up[unTracks];
949  Int_t unp(0);
950  Int_t unn(0);
951 // if(fDebug > 1) cout << " started with " << unTracks << " of tracks "<< endl;
952  for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
953  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
954  if(!track) AliFatal("Not a standard AOD");
955  if (!PhiTrack(track)) continue;
956  qx+=TMath::Cos(fHarmonic*track->Phi());
957  qy+=TMath::Sin(fHarmonic*track->Phi());
958  if (fQA) {
959  if(track->Charge() > 0) {fEventStats->Fill(1, fCentralityWeight); fPtP->Fill(track->Pt(), fCentralityWeight);}
960  if(track->Charge() < 0) {fEventStats->Fill(2, fCentralityWeight); fPtN->Fill(track->Pt(), fCentralityWeight);}
961  }
962  if (IsKaon(track)) {
963  if(fEventMixing) MixingCandidates->Add(new AliPhiMesonHelperTrack(track->Eta(), track->Phi(), track->P(),
964  track->Px(), track->Py(), track->Pz(),
965  track->Pt(), track->Charge()));
966  if (track->Charge() > 0) {
967  up[unp] = track;
968  unp++;
969  if(fQA) {fEventStats->Fill(3, fCentralityWeight);fPtKP->Fill(track->Pt(), fCentralityWeight);}
970  }
971  else if (track->Charge() < 0) {
972  un[unn] = track;
973  unn++;
974  if(fQA) {fEventStats->Fill(4, fCentralityWeight); fPtKN->Fill(track->Pt(), fCentralityWeight);}
975  }
976  }
977  }
978  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
979  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
980  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
981  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
982  PtSelector(0, up[pTracks], un[nTracks]);
983  Double_t pt = PhiPt(up[pTracks], un[nTracks]);
984  Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
985  TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
986  TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
987  TVector3 c = a + b;
988  Double_t phi = c.Phi();
989  Double_t eta = c.Eta();
990  Double_t p = TMath::Sqrt(c.Px()*c.Px()+c.Py()*c.Py()+c.Pz()*c.Pz());
991  Int_t nIDs[2];
992  nIDs[0] = up[pTracks]->GetID();
993  nIDs[1] = un[nTracks]->GetID();
994  MakeTrack(mass, pt, phi, eta, 2, nIDs, p, c.Pz());
995  }
996  }
997 // if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
998  for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
999  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
1000  if (!cand) continue;
1001 // if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
1002  for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
1003  if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
1004  for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
1005  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
1006  if (!iRP) continue;
1007  if (!iRP->InRPSelection()) continue;
1008  if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
1009 // if (fDebug > 1) printf(" was in RP set");
1010  iRP->SetForRPSelection(kFALSE);
1012  }
1013  }
1014 // if (fDebug > 1) printf("\n");
1015  }
1016  cand->SetForPOISelection(kTRUE);
1017  fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
1018  }
1019 // if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
1020  if(!fEventMixing) { // combinatorial background
1021  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
1022  for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
1023  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
1024  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
1025  PtSelector(1, up[pTracks], up[nTracks]);
1026  }
1027  }
1028  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
1029  for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
1030  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
1031  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
1032  PtSelector(2, un[nTracks], un[pTracks]);
1033  }
1034  }
1035  }
1036  if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
1037 
1038 
1039 
1040  // ugly, ugly hack but i want to pass the centrality weight to the sp task
1042  // and second guly hack, save the tpc ep
1043  fFlowEvent->SetPsi2(qx);
1044  fFlowEvent->SetPsi3(qy);
1045  PostData(1, fOutputList);
1046  PostData(2, fFlowEvent);
1047  }
1048 }
1049 //_____________________________________________________________________________
1051 {
1052  // skip the event selection for SE task (e.g. for MC productions)
1054  // exec of task se will do event selection and call UserExec
1055  else AliAnalysisTaskSE::Exec(c);
1056 }
1057 //_____________________________________________________________________________
1059 {
1060  // perform phi reconstruction with event mixing
1061 // if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
1062  AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
1063  if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
1064  if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
1065  Int_t nEvents = pool->GetCurrentNEvents();
1066 // if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
1067  for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
1068  TObjArray* mixed_candidates = pool->GetEvent(iEvent);
1069  if(!mixed_candidates) continue; // this should NEVER happen
1070  Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
1071  Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
1072 // if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
1073  AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
1074  AliPhiMesonHelperTrack* buffer_up[bufferTracks];
1075  Int_t buffer_unp(0);
1076  Int_t buffer_unn(0);
1077  AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
1078  AliPhiMesonHelperTrack* mix_up[candidates];
1079  Int_t mix_unp(0);
1080  Int_t mix_unn(0);
1081  for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
1082  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
1083  if(!track) continue;
1084  if (track->Charge() > 0) {
1085  mix_up[mix_unp] = track;
1086  mix_unp++;
1087  }
1088  else if (track->Charge() < 0 ) {
1089  mix_un[mix_unn] = track;
1090  mix_unn++;
1091  }
1092  }
1093  for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
1094  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
1095  if(!track) continue;
1096  if (track->Charge() > 0) {
1097  buffer_up[buffer_unp] = track;
1098  buffer_unp++;
1099  }
1100  else if (track->Charge() < 0 ) {
1101  buffer_un[buffer_unn] = track;
1102  buffer_unn++;
1103  }
1104  }
1105  for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
1106 // if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
1107  if(!fTypeMixing) { // mix with like-sign kaons
1108  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1109 // if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1110  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
1111  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
1112  PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
1113  }
1114  }
1115  else if(fTypeMixing) { // mix with unlike-sign kaons
1116  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1117 // if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1118  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
1119  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
1120  PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
1121  }
1122  }
1123  }
1124  for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
1125 // if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
1126  if(!fTypeMixing) { // mix with like-sign kaons
1127  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1128 // if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1129  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
1130  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
1131  PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
1132  }
1133  }
1134  else if(fTypeMixing) { // mix with unlike-sign kaons
1135  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1136 // if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1137  // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
1138  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
1139  PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
1140  }
1141  }
1142  }
1143  } // end of mixed events loop
1144  } // end of checking to see whether pool is filled correctly
1145  pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
1146 // if(fDebug > 0 ) pool->PrintInfo();
1147 }
1148 //_____________________________________________________________________________
1150 {
1151  // terminate
1152  if(fDebug > 0) cout << " *** Terminate() *** " << endl;
1153 }
1154 //______________________________________________________________________________
1156 {
1157  // Construct Flow Candidate Track from two selected candidates
1158 // if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
1159  // if requested, check rapidity at this point
1160  if(fCandidateYCut) {
1161  Double_t y = 0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
1162  if (y > fCandidateMaxY || y < fCandidateMinY) return;
1163  }
1164  Bool_t overwrite = kTRUE;
1165  AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
1166  if (!sTrack) {
1167  sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
1168  overwrite = kFALSE;
1169  }
1170  else sTrack->ClearMe();
1171  sTrack->SetMass(mass);
1172  sTrack->SetPt(pt);
1173  sTrack->SetPhi(phi);
1174  sTrack->SetEta(eta);
1175  for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
1176  sTrack->SetForPOISelection(kTRUE);
1177  sTrack->SetForRPSelection(kFALSE);
1178  if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
1179  else fCandidates->AddLast(sTrack);
1180  return;
1181 }
1182 //_____________________________________________________________________________
1184 {
1185  // Fill QA histos for MC analysis
1186  TClonesArray *arrayMC = 0;
1187  if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
1188  // fill array with mc tracks
1189  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1190  if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
1191  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
1192  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
1193  if(!track) AliFatal("Not a standard AOD");
1194  if(!PhiTrack(track) || !IsKaon(track)) { // check for kaons
1195  if(fDebug>1) cout << " Rejected track" << endl;
1196  continue;
1197  }
1198  if (fDebug>1) cout << " Received MC kaon " << endl;
1199  Double_t b[2] = { -99., -99.};
1200  Double_t bCov[3] = { -99., -99., -99.};
1201  AliAODTrack copy(*track);
1202  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
1203  // find corresponding mc particle
1204  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
1205  if (!partMC) {
1206  if(fDebug > 1) cout << "Cannot get MC particle" << endl;
1207  continue;
1208  }
1209  // Check if it is primary, secondary from material or secondary from weak decay
1210  Bool_t isPrimary = partMC->IsPhysicalPrimary();
1211  Bool_t isSecondaryMaterial = kFALSE;
1212  Bool_t isSecondaryWeak = kFALSE;
1213  if (!isPrimary) {
1214  Int_t mfl = -999, codemoth = -999;
1215  Int_t indexMoth = partMC->GetMother();
1216  if (indexMoth >= 0) { // is not fake
1217  AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
1218  codemoth = TMath::Abs(moth->GetPdgCode());
1219  mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1220  }
1221  if (mfl == 3) isSecondaryWeak = kTRUE;
1222  else isSecondaryMaterial = kTRUE;
1223  }
1224  if (isPrimary) {
1225  fDCAPrim->Fill(track->Pt(), b[0]);
1226  fDCAXYQA->Fill(b[0]);
1227  fDCAZQA->Fill(b[1]);
1228  }
1229  if (isSecondaryWeak) fDCASecondaryWeak->Fill(track->Pt(), b[0]);
1230  if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
1231  }
1232 }
1233 
1234 
1235 
1236 //_____________________________________________________________________________
1237 Double_t AliAnalysisTaskPhiFlow::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1238 {
1239 
1240  // calculate sqrt of weighted distance to other vertex
1241  if (!v0 || !v1) {
1242  printf("One of vertices is not valid\n");
1243  return 0;
1244  }
1245  static TMatrixDSym vVb(3);
1246  double dist = -1;
1247  double dx = v0->GetX()-v1->GetX();
1248  double dy = v0->GetY()-v1->GetY();
1249  double dz = v0->GetZ()-v1->GetZ();
1250  double cov0[6],cov1[6];
1251  v0->GetCovarianceMatrix(cov0);
1252  v1->GetCovarianceMatrix(cov1);
1253  //
1254  // fQxavsV0[0] fQxnmV0A
1255  // fQyavsV0[0] fQynmV0A
1256  // fQxavsV0[1] fQxnsV0A
1257  // fQyavsV0[1] fQynsV0A
1258  // fQxcvsV0[0] fQxnmV0C
1259  // fQycvsV0[0] fQynmV0C
1260  // fQxcvsV0[1] fQxnsV0C
1261  // fQycvsV0[1] fQynsV0C vVb(0,0) = cov0[0]+cov1[0];
1262  vVb(1,1) = cov0[2]+cov1[2];
1263  vVb(2,2) = cov0[5]+cov1[5];
1264  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1265  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1266  vVb.InvertFast();
1267  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1268  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1269  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1270  return dist>0 ? TMath::Sqrt(dist) : -1;
1271 
1272 }
1273 
1274 
1275 
1276 
1277 //_____________________________________________________________________________
1279 {
1280  // check for multi-vertexer pile-up
1281  //
1282  const int kMinPlpContrib = 5;
1283  const double kMaxPlpChi2 = 5.0;
1284  const double kMinWDist = 15;
1285  //
1286  const AliVVertex* vtPrm = 0;
1287  const AliVVertex* vtPlp = 0;
1288  int nPlp = 0;
1289  //
1290  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1291  vtPrm = aod->GetPrimaryVertex();
1292  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1293 
1294  //int bcPrim = vtPrm->GetBC();
1295  //
1296  for (int ipl=0;ipl<nPlp;ipl++) {
1297  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1298  //
1299  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1300  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1301  // int bcPlp = vtPlp->GetBC();
1302  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1303  //
1304  double wDst = GetWDist(vtPrm,vtPlp);
1305  if (wDst<kMinWDist) continue;
1306  //
1307  return kTRUE; // pile-up: well separated vertices
1308  }
1309  //
1310  return kFALSE;
1311  //
1312 }
1313 
1314 //_____________________________________________________________________________
1315 //---------------------------------------------------
1317 {
1318 
1319  Short_t kPID = 0;
1320 
1321  if((nSk == nSpi) && (nSk == nSp))
1322  return kPID;
1323  if((nSk < nSpi) && (nSk < nSp) && (nSk < fPIDConfig[0]))
1324  kPID = 2;
1325 
1326  if((nSpi < nSk) && (nSpi < nSp) && (nSpi < fPIDConfig[0]))
1327  kPID = 1;
1328 
1329  if((nSp < nSk) && (nSp < nSpi) && (nSp < fPIDConfig[0]))
1330  kPID = 3;
1331  return kPID;
1332 
1333 }
1334 
1335 
1336 //_____________________________________________________________________________
1338 {
1339 
1340  Bool_t kDC = kFALSE;
1341 
1342  if (nSk < fPIDConfig[0] && minNSigma != 2)
1343  kDC = kTRUE;
1344 
1345  return kDC;
1346 
1347 }
1348 
1349 
TH1F * fEta
QA plot of p_t sectrum of tracks used for event plane estimation.
Bool_t plpMV(const AliAODEvent *aod)
TH1F * fVZEROA
QA plot of eta distribution of tracks used for event plane estimation.
const Color_t cc[]
Definition: DrawKs.C:1
Bool_t PassesTPCbayesianCut(T *track) const
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)
Definition: External.C:236
TH2F * fDCAAll
QA plot TPC multiplicity (tracks used for event plane estimation)
TH1F * fPtSpectra[18]
like-sign kaon pairs
void ReconstructionWithEventMixing(TObjArray *MixingCandidates) const
Double_t PhiPt(const T *track_1, const T *track_2) const
Int_t GetNumberOfRPs() const
TH1F * fInvMNP[18]
QA histo of TOF response kaons.
Double_t mass
void AddDaughter(Int_t value)
void PlotMultiplcities(const T *event) const
void ComputeProb(const AliESDtrack *t, Float_t)
AliFlowEvent * fFlowEvent
pid response object
TCanvas * c
Definition: TestFitELoss.C:172
TH2F * fPIDk
QA histogram of TPC response of all charged particles.
TH1F * fDCAXYQA
qa dca of all charged particles
void SetMass(Double_t mass)
TH2F * fDCASecondaryWeak
dca of primaries (mc) or kaons (data)
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)
void SetReferenceMultiplicity(Int_t m)
TH1F * fPtN
QA histogram of p_t distribution of positive particles.
AliESDv0 * fV0
placeholder for the current kink
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Bool_t InRPSelection() const
AliFlowBayesianPID * fBayesianResponse
TH2F * fNOPIDTOF
QA histogram of TPC response of kaons.
AliFlowTrack * GetTrack(Int_t i)
TList * fOutputList
event pool manager
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
TH1F * fInvMNN[18]
unlike sign kaon pairs
void SetMC(Bool_t flag)
TH2F * fNOPID
QA histogram of events that do not pass centrality cut.
Bool_t fSkipEventSelection
profiles for vzero vn(minv)
void SetPsi5(Double_t gPsi5)
void SetForRPSelection(Bool_t b=kTRUE)
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
TH2F * fPIDTOF
QA histo of TOF repsonse charged particles.
static AliFlowCommonConstants * GetMaster()
AliFlowBayesianPID * fBayesianResponse
flow events (one for each inv mass band)
TProfile * fSubEventDPhiv2
dca material (mc) all (data)
AliEventPoolManager * InitializeEventMixing()
void SetPsi2(Double_t gPsi2)
Bool_t PhiTrack(T *track) const
TH2F * fMultvsCentr
QA profile global and tpc multiplicity after outlier cut.
AliPIDResponse * fPIDResponse
TH1F * fTPCM
QA plot vzeroc multiplicity (all tracks in event)
void SetPhi(Double_t phi)
short Short_t
Definition: External.C:23
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
TH2F * fDCAPrim
qa plot of dca z
void SetPt(Double_t pt)
Int_t GetIDDaughter(Int_t value) const
ClassImp(AliAnalysisTaskPhiFlow) ClassImp(AliPhiMesonHelperTrack) AliAnalysisTaskPhiFlow
Bool_t CheckCandidateEtaPtCut(const T *track1, const T *track2) const
TObjArray * fCandidates
PID response object.
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]
Input train file.
TH1F * fPtKN
QA histogram of p_t distribution of positive kaons.
TH1F * BookHistogram(const char *name)
Bool_t GetDoubleCountingK(Double_t nSk, Short_t minNSigma) const
void SetNewTrackParam(Bool_t flag=kTRUE)
TH1F * fDCAZQA
qa plot of dca xz
TH1F * InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
void SetForPOISelection(Bool_t b=kTRUE)
AliEventPoolManager * fPoolManager
AOD oject.
TH1F * fInvMPP[18]
like-sign kaon pairs
TH2F * fDCAMaterial
dca of weak (mc)
TH2F * fMultCorAfterCuts
QA histogram of p_t distribution of negative kaons.
virtual Int_t Charge() const
const char Option_t
Definition: External.C:48
Bool_t IsKaon(AliAODTrack *track) const
bool Bool_t
Definition: External.C:53
virtual void Terminate(Option_t *)
TH1F * fVZEROC
QA plot vzeroa multiplicity (all tracks in event)
TH1F * BookPtHistogram(const char *name)
Short_t FindMinNSigma(Double_t nSpi, Double_t nSk, Double_t nSp) const
TH2F * BookPIDHistogram(const char *name, Bool_t TPC)
TH1F * fPtKP
QA histogram of p_t distribution of negative particles.
TProfile * fV0Data[18][2]
subevent resolution info for v2
virtual void UserExec(Option_t *option)
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
Definition: External.C:196
TH1F * fPt
QA plot of azimuthal distribution of tracks used for event plane estimation.
void PtSelector(Int_t _track_type, const T *track_1, const T *track_2) const
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