AliPhysics  01c3411 (01c3411)
 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 
61 ClassImp(AliAnalysisTaskPhiFlow)
62 ClassImp(AliPhiMesonHelperTrack)
63 
64 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
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), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0), fPileUp(kTRUE)
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 //_____________________________________________________________________________
81 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
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), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0), fPileUp(kTRUE)
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 //_____________________________________________________________________________
102 AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
103 {
104  // Destructor
105  if (fNullCuts) delete fNullCuts;
106  if (fOutputList) delete fOutputList;
107  if (fCandidates) delete fCandidates;
108  if (fFlowEvent) delete fFlowEvent;
109  if (fBayesianResponse) delete fBayesianResponse;
110  if (fEventMixing) delete fPoolManager;
111  if (fPIDCombined) delete fPIDCombined;
112  if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskPhiFlow === " << endl;
113 }
114 //_____________________________________________________________________________
115 TH1F* AliAnalysisTaskPhiFlow::BookHistogram(const char* name)
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 //_____________________________________________________________________________
128 TH2F* AliAnalysisTaskPhiFlow::BookPIDHistogram(const char* name, Bool_t TPC)
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 //_____________________________________________________________________________
146 TH1F* AliAnalysisTaskPhiFlow::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
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 //_____________________________________________________________________________
156 TH1F* AliAnalysisTaskPhiFlow::BookPtHistogram(const char* name)
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 //_____________________________________________________________________________
168 void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
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);
179  fOutputList->Add(fCentralityPass);
180  if(fQA) {
181  fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
182  fOutputList->Add(fCentralityNoPass);
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  if(fCentralityCut2010 || fCentralityCut2011) {
216  fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
217  fOutputList->Add(fMultCorAfterCuts);
218  fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
219  fOutputList->Add(fMultvsCentr);
220  }
221  }
222  if(fIsMC || fQA) {
223  fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
224  fOutputList->Add(fDCAAll);
225  fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
226  fOutputList->Add(fDCAPrim);
227  fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
228  fOutputList->Add(fDCASecondaryWeak);
229  fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
230  fOutputList->Add(fDCAMaterial);
231  }
232  if(fV0) {
233  fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 5, 0, 5);
234  fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<#Psi_{a} - #Psi_{b}>");
235  fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<#Psi_{a} - #Psi_{c}>");
236  fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<#Psi_{b} - #Psi_{c}>");
237  fSubEventDPhiv2->GetXaxis()->SetBinLabel(4, "#sqrt{#frac{<#Psi_{a} - #Psi_{b}><#Psi_{a} - #Psi_{c}>}{<#Psi_{b} - #Psi_{c}>}}");
238  fSubEventDPhiv2->GetXaxis()->SetBinLabel(5, "#sqrt{#frac{<#Psi_{a} - #Psi_{c}><#Psi_{b} - #Psi_{c}>}{<#Psi_{a} - #Psi_{b}>}}");
239  fOutputList->Add(fSubEventDPhiv2);
240  const char* V0[] = {"V0A", "V0C"};
241  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
242  for(Int_t iV0(0); iV0 < 2; iV0++) {
243  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);
244  fOutputList->Add(fV0Data[ptbin][iV0]);
245  }
246  }
247 }
248 //_____________________________________________________________________________
249 void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
250 {
251  // Create user defined output objects
252  if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
253  fNullCuts = new AliFlowTrackCuts("null_cuts");
254  fBayesianResponse = new AliFlowBayesianPID();
255  // combined pid
256  fPIDCombined = new AliPIDCombined;
257  fPIDCombined->SetDefaultTPCPriors();
258  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
259 
260  // flag to mc
261  if(fSkipEventSelection || fIsMC) fBayesianResponse->SetMC(kTRUE);
262  Double_t t(0);
263  for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
264  if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
265  if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
266  AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
267  cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
268  cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
269  cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
270  cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
271  cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
272  cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
273  fBayesianResponse->SetNewTrackParam();
274  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
275  if (man) {
276  AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
277  if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
278  }
279  // Create all output objects and store them to a list
280  fOutputList = new TList();
281  fOutputList->SetOwner(kTRUE);
282  // Create and post the Phi identification output objects
283  AddPhiIdentificationOutputObjects();
284  PostData(1, fOutputList);
285  // create candidate array
286  fCandidates = new TObjArray(1000);
287  fCandidates->SetOwner(kTRUE);
288  // create and post flowevent
289  fFlowEvent = new AliFlowEvent(10000);
290  PostData(2, fFlowEvent);
291  if(fEventMixing) fPoolManager = InitializeEventMixing();
292 }
293 //_____________________________________________________________________________
294 AliEventPoolManager* AliAnalysisTaskPhiFlow::InitializeEventMixing()
295 {
296  // initialize event mixing
297  if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
298  Int_t _c(0), _v(0);
299  for(Int_t i(0); i < 19; i++) {
300  if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
301  else _c = 19;
302  }
303  for(Int_t i(0); i < 19; i++) {
304  if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
305  else _v = 19;
306  }
307  if(fDebug > 0 ) cout << Form(" --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
308  Double_t centralityBins[_c];
309  Double_t vertexBins[_v];
310  for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
311  for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
312  return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
313 }
314 //_____________________________________________________________________________
315 template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* track1, const T* track2) const
316 {
317  // Return the invariant mass of two tracks, assuming both tracks are kaons
318  if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
319  if ((!track2) || (!track1)) return 0.;
320  Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
321  Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
322  Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
323  Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
324  Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
325  Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
326  Double_t es = TMath::Power((e1 + e2), 2);
327  if ((es - (pxs + pys + pzs)) < 0) return 0.;
328  return TMath::Sqrt((es - (pxs + pys + pzs)));
329 }
330 //_____________________________________________________________________________
331 /*
332 template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
333 {
334  // Calculate the delta dip angle between two particles
335  if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
336  if (track1->P()*track2->P() == 0) return 999;
337  return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
338 }
339 //_____________________________________________________________________________
340 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T* track1, const T* track2) const
341 {
342  // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
343  if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
344  if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
345  return kTRUE;
346 }
347 */
348 //_____________________________________________________________________________
349 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
350 {
351  // Check if pair passes eta and pt cut
352  if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
353  if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
354  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
355  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
356  TVector3 c = a + b;
357  if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
358  return kTRUE;
359 }
360 //_____________________________________________________________________________
361 template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
362 {
363  // Impose event cuts
364  if(fDebug > 0) cout << " *** EventCut() *** " << endl;
365  if (!event) return kFALSE;
366 // if (fSkipEventSelection) return kTRUE;
367  if (!CheckVertex(event)) return kFALSE;
368  if (!CheckCentrality(event)) return kFALSE;
369  if(fQA) PlotMultiplcities(event);
370  return kTRUE;
371 }
372 //_____________________________________________________________________________
373 template <typename T> void AliAnalysisTaskPhiFlow::PlotMultiplcities(const T* event) const
374 {
375  // QA multiplicity plots
376  if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
377  fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
378  fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
379  fTPCM->Fill(event->GetNumberOfTracks());
380 }
381 //_____________________________________________________________________________
382 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckVertex(const T* event)
383 {
384  // Check if event vertex is within given range
385  if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
386  if (!event->GetPrimaryVertex()) return 0x0;
387  fVertex = event->GetPrimaryVertex()->GetZ();
388  if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
389  return kTRUE;
390 }
391 //_____________________________________________________________________________
392 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
393 {
394  // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
395  if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
396  if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
397 
398  // check if the AliMultSelection object is present. If so, we should invoke the
399  // new centrality framework
400 
401  AliMultSelection *multSelection = 0x0;
402  multSelection = static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
403  if(multSelection) {
404  fCentrality = multSelection->GetMultiplicityPercentile("V0M");
405  if(fCentrality > fCentralityMin && fCentrality < fCentralityMax) {
406  fCentralityPass->Fill(fCentrality);
407  return kTRUE;
408  } else {
409  fCentralityNoPass->Fill(fCentrality);
410  return kFALSE;
411  }
412  }
413 
414  else fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
415  Double_t cenB(-999);
416  // if a second centrality estimator is requited, set it
417  (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
418  if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
419  if(fQA) fCentralityNoPass->Fill(fCentrality) ;
420  return kFALSE;
421  }
422  const Int_t nGoodTracks = event->GetNumberOfTracks();
423  if(fCentralityCut2010) { // cut on outliers
424  Float_t multTPC(0.); // tpc mult estimate
425  Float_t multGlob(0.); // global multiplicity
426  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
427  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
428  if(!trackAOD) AliFatal("Not a standard AOD");
429  if (!trackAOD) continue;
430  if (!(trackAOD->TestFilterBit(1))) continue;
431  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;
432  multTPC++;
433  }
434  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
435  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
436  if(!trackAOD) AliFatal("Not a standard AOD");
437  if (!trackAOD) continue;
438  if (!(trackAOD->TestFilterBit(16))) continue;
439  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;
440  Double_t b[2] = {-99., -99.};
441  Double_t bCov[3] = {-99., -99., -99.};
442  AliAODTrack copy(*trackAOD);
443  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
444  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
445  multGlob++;
446  } //track loop
447  // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
448  if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
449  if(fQA) {
450  fMultCorAfterCuts->Fill(multGlob, multTPC);
451  fMultvsCentr->Fill(fCentrality, multTPC);
452  }
453  }
454 
455  if(fCentralityCut2011) { // 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 > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
481  if(fQA) {
482  fMultCorAfterCuts->Fill(multGlob, multTPC);
483  fMultvsCentr->Fill(fCentrality, multTPC);
484  }
485  }
486 
487  fCentralityPass->Fill(fCentrality);
488  return kTRUE;
489 }
490 //_____________________________________________________________________________
491 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
492 {
493  // Initialize the Bayesian PID object for AOD
494  if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
495  fBayesianResponse->SetDetResponse(event, fCentrality);
496 }
497 //_____________________________________________________________________________
498 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PassesTPCbayesianCut(T* track) const
499 {
500  // Check if the particle passes the TPC TOF bayesian cut.
501  if(fDebug > 1) cout << " *** PassesTPCbayesianCut() *** " << endl;
502  fBayesianResponse->ComputeProb(track);
503  if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
504  Float_t *probabilities = fBayesianResponse->GetProb();
505  if (probabilities[3] > fPIDConfig[6]) {
506  if(fQA) {fPhi->Fill(track->Phi()); fPt->Fill(track->Pt()); fEta->Fill(track->Eta());}
507  return kTRUE;
508  }
509  return kFALSE;
510 }
511 //_____________________________________________________________________________
512 Bool_t AliAnalysisTaskPhiFlow::PassesDCACut(AliAODTrack* track) const
513 {
514  // check if track passes dca cut according to dca routine
515  // setup the routine as follows:
516  // fDCAConfig[0] < -1 no pt dependence
517  // fDCAConfig[0] = 0 do nothing
518  // fDCAConfig[0] > 1 pt dependent dca cut
519  if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
520  if(fIsMC) return kTRUE;
521  AliAODTrack copy(*track);
522  if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
523  if(fQA) {
524  Double_t b[2] = { -99., -99.};
525  Double_t bCov[3] = { -99., -99., -99.};
526  if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
527  fDCAXYQA->Fill(b[0]);
528  fDCAZQA->Fill(b[1]);
529  fDCAPrim->Fill(track->Pt(), b[0]);
530  }
531  }
532  return kTRUE;
533  }
534  Double_t b[2] = { -99., -99.};
535  Double_t bCov[3] = { -99., -99., -99.};
536  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
537  if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
538  if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
539  if(fDCAConfig[0] > .9) {
540  if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
541  Double_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
542  if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
543  if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
544  }
545  if(fQA) {
546  fDCAXYQA->Fill(b[0]);
547  fDCAZQA->Fill(b[1]);
548  fDCAPrim->Fill(track->Pt(), b[0]);
549  }
550  return kTRUE;
551 }
552 //_____________________________________________________________________________
553 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
554 {
555  // Kaon identification routine, based on multiple detectors and approaches
556  if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
557 
558 /* if(fUsePidResponse) {
559  Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
560  fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
561  if(prob[3] > fPIDConfig[6]) return kTRUE;
562  }
563  if(!PassesDCACut(track)) return kFALSE;
564  if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
565  if(track->Pt() < fPIDConfig[1]) {
566  if(fDebug>1) cout << " ITS received track with p_t " << track->Pt() << endl;
567  // if tpc control is disabled, pure its pid
568  if(fPIDConfig[2] < 0.) {
569  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) return kTRUE;
570  return kFALSE;
571  }
572  // else, switch to ITS pid with TPC rejection of protons and pions
573  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton)) < fPIDConfig[3]) return kFALSE;
574  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion)) < fPIDConfig[3]) return kFALSE;
575  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kKaon)) < fPIDConfig[0]) {
576  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
577  return kTRUE;
578  }
579  return kFALSE;
580  }
581  if((track->Pt() > fPIDConfig[1]) && (track->Pt() < fPIDConfig[4])) {
582  if(fDebug>1) cout << " TPC received track with p_t " << track->Pt() << endl;
583  // if its control is disabled, pure tpc pid
584  if(fPIDConfig[5] < 0.) {
585  if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) return kTRUE;
586  return kFALSE;
587  }
588  // else, switch to TPC pid with ITS rejection of protons and pions
589  if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kProton)) < fPIDConfig[0]) return kFALSE;
590  else if (TMath::Abs(fPIDResponse->NumberOfSigmasITS(track, AliPID::kPion)) < fPIDConfig[0]) return kFALSE;
591  else if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon)) < fPIDConfig[3]) {
592  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal()); fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
593  return kTRUE;
594  }
595  return kFALSE;
596  }
597  if(fDebug>1) cout << " Bayesian method received track with p_t " << track->Pt() << endl;
598  // switch to bayesian PID
599  if (PassesTPCbayesianCut(track)) {
600  if(fQA) {fPIDk->Fill(track->P(), track->GetTPCsignal());fPIDTOF->Fill(track->P(), track->GetTOFsignal());}
601  return kTRUE;
602  }*/
603  if(fQA) {
604  Float_t length = track->GetIntegratedLength();
605  Float_t time = track->GetTOFsignal();
606  Double_t beta = -.05;
607 
608  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
609 
610  fNOPID->Fill(track->P(), track->GetTPCsignal());
611  fNOPIDTOF->Fill(track->P(), beta);
612  }
613 
614 
615 
616  //kaon stuff
617  Double_t nSigKTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kKaon);
618  Double_t nSigKTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kKaon);
619 
620  Double_t nSigmaK = 999;
621 
622  if(track->Pt() > .4) {
623  nSigmaK = TMath::Sqrt(nSigKTPC*nSigKTPC+nSigKTOF*nSigKTOF);
624  } else {
625  if(track->GetTPCsignal() > 110) nSigmaK = TMath::Abs(nSigKTPC);
626  }
627  // pion
628  Double_t nSigPiTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
629  Double_t nSigPiTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kPion);
630 
631  Double_t nSigmaPi = 999;
632 
633  if(track->Pt() > .4) {
634  nSigmaPi = TMath::Sqrt(nSigPiTPC*nSigPiTPC+nSigPiTOF*nSigPiTOF);
635  } else {
636  if(track->GetTPCsignal() < 60) nSigmaPi = TMath::Abs(nSigPiTPC);
637  }
638 
639  //proton
640  Double_t nSigPTPC = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton);
641  Double_t nSigPTOF = fPIDResponse->NumberOfSigmasTOF(track, AliPID::kProton);
642 
643  Double_t nSigmaP = 999;
644 
645  if(track->Pt() > .4) {
646  nSigmaP = TMath::Sqrt(nSigPTPC*nSigPTPC+nSigPTOF*nSigPTOF);
647  } else {
648  if(track->GetTPCsignal() > 110) nSigmaP = TMath::Abs(nSigPTPC);
649  }
650 
651  Short_t minSigma = FindMinNSigma(nSigmaPi, nSigmaK, nSigmaP);
652 
653 
654 
655  if(minSigma == 0) return kFALSE;
656  if((nSigmaK == nSigmaPi) && ( nSigmaK == nSigmaP)) return kFALSE;
657 
658  if(minSigma == 2 &&!GetDoubleCountingK(nSigmaK, minSigma)) {
659 
660  if(fQA) {
661  Float_t length = track->GetIntegratedLength();
662  Float_t time = track->GetTOFsignal();
663  Double_t beta = -.05;
664 
665  if((length > 0) && (time > 0)) beta = length / 2.99792458e-2 / time;
666 
667  // additional selection to bruteforce throw away bad betas
668  if(beta < 0.4) return kFALSE;
669 
670  fPIDk->Fill(track->P(), track->GetTPCsignal());
671  if(track->Pt() > .4) fPIDTOF->Fill(track->P(), beta);
672  }
673 
674 
675  return kTRUE;
676  }
677 
678  return kFALSE;
679 }
680 //_____________________________________________________________________________
681 template <typename T> Double_t AliAnalysisTaskPhiFlow::PhiPt(const T* track1, const T* track2) const
682 {
683  // return p_t of track pair
684  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
685  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
686  TVector3 c = a + b;
687  return c.Pt();
688 }
689 //_____________________________________________________________________________
690 template <typename T> void AliAnalysisTaskPhiFlow::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
691 {
692  // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
693  Double_t pt = PhiPt(track1, track2);
694  if (tracktype == 0) {
695  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
696  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
697  fInvMNP[ptbin]->Fill(InvariantMass(track1, track2));
698  if(fQA) fPtSpectra[ptbin]->Fill(pt);
699  }
700  }
701  }
702  if (tracktype == 1) {
703  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
704  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
705  fInvMPP[ptbin]->Fill(InvariantMass(track1, track2));
706  }
707  }
708  }
709  if (tracktype == 2) {
710  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
711  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
712  fInvMNN[ptbin]->Fill(InvariantMass(track1, track2));
713  }
714  }
715  }
716 }
717 //_____________________________________________________________________________
718 template <typename T> Bool_t AliAnalysisTaskPhiFlow::PhiTrack(T* track) const
719 {
720  // check if track meets quality cuts
721  if(!track) return kFALSE;
722  return fPOICuts->IsSelected(track);
723 }
724 //_____________________________________________________________________________
725 template <typename T> void AliAnalysisTaskPhiFlow::SetNullCuts(T* event)
726 {
727  // Set null cuts
728  if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
729  fCutsRP->SetEvent(event, MCEvent());
730  fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
731  fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
732  fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
733  fNullCuts->SetEvent(event, MCEvent());
734 }
735 //_____________________________________________________________________________
736 void AliAnalysisTaskPhiFlow::PrepareFlowEvent(Int_t iMulti)
737 {
738  // Prepare flow events
739  if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
740  fFlowEvent->ClearFast();
741  fFlowEvent->Fill(fCutsRP, fNullCuts);
742  fFlowEvent->SetReferenceMultiplicity(iMulti);
743  fFlowEvent->DefineDeadZone(0, 0, 0, 0);
744 }
745 //_____________________________________________________________________________
746 void AliAnalysisTaskPhiFlow::VZEROSubEventAnalysis()
747 {
748  // vzero event plane analysis using three subevents
749  if(fDebug > 0) cout << " ** VZEROSubEventAnalysis() *** " << endl;
750  if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
751  if(fDebug > 0 ) cout << "Fatal error: unable to retrieve VZERO task output !!! Exiting ..." << endl;
752  return;
753  }
754  // retrieve data
755  Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
756  for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0) {
757  if(fDebug > 0) cout << " Warning: I've encountered 0 in one of the symmetry panes (TPC, VZERO) - skipping event !" << endl;
758  return;
759  }
760  // save info for resolution calculations
761  Float_t qaqb = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1])); // vzeroa - tpc
762  Float_t qaqc = TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2])); // vzeroa - vzeroc
763  Float_t qbqc = TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2])); // tpc - vzeroc
764  fSubEventDPhiv2->Fill(0.5, qaqb);
765  fSubEventDPhiv2->Fill(1.5, qaqc);
766  fSubEventDPhiv2->Fill(2.5, qbqc);
767  Float_t minv[31];
768  Float_t _dphi[30][fNPtBins][2]; //minv, pt, v0a-c
769  Int_t occurence[30][fNPtBins]; //minv, pt
770  for(Int_t mb(0); mb < 31; mb++) minv[mb] = 0.99 + mb * 0.0034;
771  for(Int_t i(0); i < 30; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
772  _dphi[i][j][k] = 0;
773  if(k==0) occurence[i][j] = 0;
774  }
775  for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign kaon pairs
776  AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
777  if(!track) {
778  if(fDebug > 1) cout << " --> dynamic_cast returned null-pointer, skipping candidate <-- " << endl;
779  continue;
780  }
781  for(Int_t mb(0); mb < 30; mb++) { // loop over massbands
782  if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
783  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
784  if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
785  _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
786  _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
787  occurence[mb][ptbin]+=1;
788  }
789  }
790  for(Int_t mb(0); mb < 30; mb++) // write vn values to tprofiles
791  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
792  if(occurence[mb][ptbin]==0) continue;
793  fV0Data[ptbin][0]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
794  fV0Data[ptbin][1]->Fill(mb*0.0034+0.99+0.0017, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
795  }
796  }
797 }
798 //_____________________________________________________________________________
799 void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
800 {
801  // UserExec: called for each event. Commented where necessary
802  if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
803  TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
804  if(fEventMixing) {
805  MixingCandidates = new TObjArray();
806  MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
807  }
808  if (!fPIDResponse) {
809  if(fDebug > 0 ) cout << " Could not get PID response " << endl;
810  return;
811  }
812  fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
813  if (fAOD) {
814  if (!EventCut(fAOD)) return; // check for event cuts
815 
816 
817  // add extra pileup cuts for high intensity runs
818  // courtesy of alex dobrin
819  if (fPileUp){
820 
821  if (plpMV(fAOD))
822  return;
823 
824  Short_t isPileup = fAOD->IsPileupFromSPD(3);
825  if (isPileup != 0)
826  return;
827 
828  if (((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0)
829  return;
830 
831  /*
832  Int_t bc2 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
833  if (bc2 != 0)
834  return;
835 
836  Int_t bc1 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
837  if (bc1 != 0)
838  return;
839  */
840 
841  // add vertexer selection
842  AliAODVertex* vtTrc = fAOD->GetPrimaryVertex();
843  AliAODVertex* vtSPD = fAOD->GetPrimaryVertexSPD();
844  if (vtTrc->GetNContributors()<2 || vtSPD->GetNContributors()<1) return; // one of vertices is missing
845  double covTrc[6],covSPD[6];
846  vtTrc->GetCovarianceMatrix(covTrc);
847  vtSPD->GetCovarianceMatrix(covSPD);
848  double dz = vtTrc->GetZ()-vtSPD->GetZ();
849  double errTot = TMath::Sqrt(covTrc[5]+covSPD[5]);
850  double errTrc = TMath::Sqrt(covTrc[5]);
851  double nsigTot = TMath::Abs(dz)/errTot, nsigTrc = TMath::Abs(dz)/errTrc;
852  if (TMath::Abs(dz)>0.2 || nsigTot>10 || nsigTrc>20) return; // bad vertexing
853 
854 
855 
856  }
857 
858  //new function for 2015 to remove incomplete events
859  if (fAOD->IsIncompleteDAQ())
860  return;
861 
862 
863 // InitializeBayesianPID(fAOD); // init event objects
864  SetNullCuts(fAOD);
865  PrepareFlowEvent(fAOD->GetNumberOfTracks());
866  fCandidates->SetLast(-1);
867  if(fIsMC) IsMC(); // launch mc stuff
868  if(fQA) fEventStats->Fill(0);
869  Int_t unTracks = fAOD->GetNumberOfTracks();
870  AliAODTrack* un[unTracks];
871  AliAODTrack* up[unTracks];
872  Int_t unp(0);
873  Int_t unn(0);
874  if(fDebug > 1) cout << " started with " << unTracks << " of tracks "<< endl;
875  for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
876  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
877  if(!track) AliFatal("Not a standard AOD");
878  if (!PhiTrack(track)) continue;
879  if (fQA) {
880  if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
881  if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
882  }
883  if (IsKaon(track)) {
884  if(fEventMixing) MixingCandidates->Add(new AliPhiMesonHelperTrack(track->Eta(), track->Phi(), track->P(),
885  track->Px(), track->Py(), track->Pz(),
886  track->Pt(), track->Charge()));
887  if (track->Charge() > 0) {
888  up[unp] = track;
889  unp++;
890  if(fQA) {fEventStats->Fill(3);fPtKP->Fill(track->Pt());}
891  }
892  else if (track->Charge() < 0) {
893  un[unn] = track;
894  unn++;
895  if(fQA) {fEventStats->Fill(4); fPtKN->Fill(track->Pt());}
896  }
897  }
898  }
899  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
900  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
901 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
902  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
903  PtSelector(0, up[pTracks], un[nTracks]);
904  Double_t pt = PhiPt(up[pTracks], un[nTracks]);
905  Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
906  TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
907  TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
908  TVector3 c = a + b;
909  Double_t phi = c.Phi();
910  Double_t eta = c.Eta();
911  Double_t p = TMath::Sqrt(c.Px()*c.Px()+c.Py()*c.Py()+c.Pz()*c.Pz());
912  Int_t nIDs[2];
913  nIDs[0] = up[pTracks]->GetID();
914  nIDs[1] = un[nTracks]->GetID();
915  MakeTrack(mass, pt, phi, eta, 2, nIDs, p, c.Pz());
916  }
917  }
918  if(fV0) VZEROSubEventAnalysis();
919  if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
920  for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
921  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
922  if (!cand) continue;
923  if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
924  for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
925  if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
926  for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
927  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
928  if (!iRP) continue;
929  if (!iRP->InRPSelection()) continue;
930  if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
931  if (fDebug > 1) printf(" was in RP set");
932  iRP->SetForRPSelection(kFALSE);
933  fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
934  }
935  }
936  if (fDebug > 1) printf("\n");
937  }
938  cand->SetForPOISelection(kTRUE);
939  fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
940  }
941  if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
942  if(!fEventMixing) { // combinatorial background
943  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
944  for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
945 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
946  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
947  PtSelector(1, up[pTracks], up[nTracks]);
948  }
949  }
950  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
951  for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
952 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
953  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
954  PtSelector(2, un[nTracks], un[pTracks]);
955  }
956  }
957  }
958  if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
959  PostData(1, fOutputList);
960  PostData(2, fFlowEvent);
961  }
962 }
963 //_____________________________________________________________________________
964 void AliAnalysisTaskPhiFlow::Exec(Option_t* c)
965 {
966  // skip the event selection for SE task (e.g. for MC productions)
967  if(fSkipEventSelection) AliAnalysisTaskPhiFlow::UserExec(c);
968  // exec of task se will do event selection and call UserExec
969  else AliAnalysisTaskSE::Exec(c);
970 }
971 //_____________________________________________________________________________
972 void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
973 {
974  // perform phi reconstruction with event mixing
975  if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
976  AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
977  if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
978  if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
979  Int_t nEvents = pool->GetCurrentNEvents();
980  if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
981  for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
982  TObjArray* mixed_candidates = pool->GetEvent(iEvent);
983  if(!mixed_candidates) continue; // this should NEVER happen
984  Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
985  Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
986  if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
987  AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
988  AliPhiMesonHelperTrack* buffer_up[bufferTracks];
989  Int_t buffer_unp(0);
990  Int_t buffer_unn(0);
991  AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
992  AliPhiMesonHelperTrack* mix_up[candidates];
993  Int_t mix_unp(0);
994  Int_t mix_unn(0);
995  for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
996  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
997  if(!track) continue;
998  if (track->Charge() > 0) {
999  mix_up[mix_unp] = track;
1000  mix_unp++;
1001  }
1002  else if (track->Charge() < 0 ) {
1003  mix_un[mix_unn] = track;
1004  mix_unn++;
1005  }
1006  }
1007  for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
1008  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
1009  if(!track) continue;
1010  if (track->Charge() > 0) {
1011  buffer_up[buffer_unp] = track;
1012  buffer_unp++;
1013  }
1014  else if (track->Charge() < 0 ) {
1015  buffer_un[buffer_unn] = track;
1016  buffer_unn++;
1017  }
1018  }
1019  for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
1020  if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
1021  if(!fTypeMixing) { // mix with like-sign kaons
1022  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1023  if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1024 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
1025  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
1026  PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
1027  }
1028  }
1029  else if(fTypeMixing) { // mix with unlike-sign kaons
1030  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1031  if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1032 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
1033  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
1034  PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
1035  }
1036  }
1037  }
1038  for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
1039  if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
1040  if(!fTypeMixing) { // mix with like-sign kaons
1041  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
1042  if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
1043 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
1044  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
1045  PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
1046  }
1047  }
1048  else if(fTypeMixing) { // mix with unlike-sign kaons
1049  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
1050  if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
1051 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
1052  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
1053  PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
1054  }
1055  }
1056  }
1057  } // end of mixed events loop
1058  } // end of checking to see whether pool is filled correctly
1059  pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
1060  if(fDebug > 0 ) pool->PrintInfo();
1061 }
1062 //_____________________________________________________________________________
1063 void AliAnalysisTaskPhiFlow::Terminate(Option_t *)
1064 {
1065  // terminate
1066  if(fDebug > 0) cout << " *** Terminate() *** " << endl;
1067 }
1068 //______________________________________________________________________________
1069 void AliAnalysisTaskPhiFlow::MakeTrack(Double_t mass, Double_t pt, Double_t phi, Double_t eta, Int_t nDau, Int_t iID[], Double_t p, Double_t pz) const
1070 {
1071  // Construct Flow Candidate Track from two selected candidates
1072  if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
1073  // if requested, check rapidity at this point
1074  if(fCandidateYCut) {
1075  Double_t y = 0.5*TMath::Log((TMath::Sqrt(mass*mass+p*p)+pz)/(TMath::Sqrt(mass*mass+p*p)-pz));
1076  if (y > fCandidateMaxY || y < fCandidateMinY) return;
1077  }
1078  Bool_t overwrite = kTRUE;
1079  AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
1080  if (!sTrack) {
1081  sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
1082  overwrite = kFALSE;
1083  }
1084  else sTrack->ClearMe();
1085  sTrack->SetMass(mass);
1086  sTrack->SetPt(pt);
1087  sTrack->SetPhi(phi);
1088  sTrack->SetEta(eta);
1089  for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
1090  sTrack->SetForPOISelection(kTRUE);
1091  sTrack->SetForRPSelection(kFALSE);
1092  if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
1093  else fCandidates->AddLast(sTrack);
1094  return;
1095 }
1096 //_____________________________________________________________________________
1097 void AliAnalysisTaskPhiFlow::IsMC()
1098 {
1099  // Fill QA histos for MC analysis
1100  TClonesArray *arrayMC = 0;
1101  if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
1102  // fill array with mc tracks
1103  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
1104  if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
1105  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
1106  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
1107  if(!track) AliFatal("Not a standard AOD");
1108  if(!PhiTrack(track) || !IsKaon(track)) { // check for kaons
1109  if(fDebug>1) cout << " Rejected track" << endl;
1110  continue;
1111  }
1112  if (fDebug>1) cout << " Received MC kaon " << endl;
1113  Double_t b[2] = { -99., -99.};
1114  Double_t bCov[3] = { -99., -99., -99.};
1115  AliAODTrack copy(*track);
1116  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
1117  // find corresponding mc particle
1118  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
1119  if (!partMC) {
1120  if(fDebug > 1) cout << "Cannot get MC particle" << endl;
1121  continue;
1122  }
1123  // Check if it is primary, secondary from material or secondary from weak decay
1124  Bool_t isPrimary = partMC->IsPhysicalPrimary();
1125  Bool_t isSecondaryMaterial = kFALSE;
1126  Bool_t isSecondaryWeak = kFALSE;
1127  if (!isPrimary) {
1128  Int_t mfl = -999, codemoth = -999;
1129  Int_t indexMoth = partMC->GetMother();
1130  if (indexMoth >= 0) { // is not fake
1131  AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
1132  codemoth = TMath::Abs(moth->GetPdgCode());
1133  mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
1134  }
1135  if (mfl == 3) isSecondaryWeak = kTRUE;
1136  else isSecondaryMaterial = kTRUE;
1137  }
1138  if (isPrimary) {
1139  fDCAPrim->Fill(track->Pt(), b[0]);
1140  fDCAXYQA->Fill(b[0]);
1141  fDCAZQA->Fill(b[1]);
1142  }
1143  if (isSecondaryWeak) fDCASecondaryWeak->Fill(track->Pt(), b[0]);
1144  if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
1145  }
1146 }
1147 
1148 
1149 
1150 //_____________________________________________________________________________
1151 Double_t AliAnalysisTaskPhiFlow::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
1152 {
1153 
1154  // calculate sqrt of weighted distance to other vertex
1155  if (!v0 || !v1) {
1156  printf("One of vertices is not valid\n");
1157  return 0;
1158  }
1159  static TMatrixDSym vVb(3);
1160  double dist = -1;
1161  double dx = v0->GetX()-v1->GetX();
1162  double dy = v0->GetY()-v1->GetY();
1163  double dz = v0->GetZ()-v1->GetZ();
1164  double cov0[6],cov1[6];
1165  v0->GetCovarianceMatrix(cov0);
1166  v1->GetCovarianceMatrix(cov1);
1167  //
1168  // fQxavsV0[0] fQxnmV0A
1169  // fQyavsV0[0] fQynmV0A
1170  // fQxavsV0[1] fQxnsV0A
1171  // fQyavsV0[1] fQynsV0A
1172  // fQxcvsV0[0] fQxnmV0C
1173  // fQycvsV0[0] fQynmV0C
1174  // fQxcvsV0[1] fQxnsV0C
1175  // fQycvsV0[1] fQynsV0C vVb(0,0) = cov0[0]+cov1[0];
1176  vVb(1,1) = cov0[2]+cov1[2];
1177  vVb(2,2) = cov0[5]+cov1[5];
1178  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1179  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1180  vVb.InvertFast();
1181  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1182  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1183  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1184  return dist>0 ? TMath::Sqrt(dist) : -1;
1185 
1186 }
1187 
1188 
1189 
1190 
1191 //_____________________________________________________________________________
1192 Bool_t AliAnalysisTaskPhiFlow::plpMV(const AliAODEvent* aod)
1193 {
1194  // check for multi-vertexer pile-up
1195  //
1196  const int kMinPlpContrib = 5;
1197  const double kMaxPlpChi2 = 5.0;
1198  const double kMinWDist = 15;
1199  //
1200  const AliVVertex* vtPrm = 0;
1201  const AliVVertex* vtPlp = 0;
1202  int nPlp = 0;
1203  //
1204  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1205  vtPrm = aod->GetPrimaryVertex();
1206  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1207 
1208  //int bcPrim = vtPrm->GetBC();
1209  //
1210  for (int ipl=0;ipl<nPlp;ipl++) {
1211  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1212  //
1213  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1214  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1215  // int bcPlp = vtPlp->GetBC();
1216  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1217  //
1218  double wDst = GetWDist(vtPrm,vtPlp);
1219  if (wDst<kMinWDist) continue;
1220  //
1221  return kTRUE; // pile-up: well separated vertices
1222  }
1223  //
1224  return kFALSE;
1225  //
1226 }
1227 
1228 //_____________________________________________________________________________
1229 //---------------------------------------------------
1230 Short_t AliAnalysisTaskPhiFlow::FindMinNSigma(Double_t nSpi, Double_t nSk, Double_t nSp) const
1231 {
1232 
1233  Short_t kPID = 0;
1234 
1235  if((nSk == nSpi) && (nSk == nSp))
1236  return kPID;
1237  if((nSk < nSpi) && (nSk < nSp) && (nSk < fPIDConfig[0]))
1238  kPID = 2;
1239 
1240  if((nSpi < nSk) && (nSpi < nSp) && (nSpi < fPIDConfig[0]))
1241  kPID = 1;
1242 
1243  if((nSp < nSk) && (nSp < nSpi) && (nSp < fPIDConfig[0]))
1244  kPID = 3;
1245  return kPID;
1246 
1247 }
1248 
1249 
1250 //_____________________________________________________________________________
1251 Bool_t AliAnalysisTaskPhiFlow::GetDoubleCountingK(Double_t nSk, Short_t minNSigma) const
1252 {
1253 
1254  Bool_t kDC = kFALSE;
1255 
1256  if (nSk < fPIDConfig[0] && minNSigma != 2)
1257  kDC = kTRUE;
1258 
1259  return kDC;
1260 
1261 }
1262 
AliFlowEvent.h
AliAnalysisTaskPhiFlow::plpMV
Bool_t plpMV(const AliAODEvent *aod)
Definition: AliAnalysisTaskPhiFlow.cxx:1192
cc
const Color_t cc[]
Definition: DrawKs.C:1
AliAnalysisTaskPhiFlow::PassesTPCbayesianCut
Bool_t PassesTPCbayesianCut(T *track) const
Definition: AliAnalysisTaskPhiFlow.cxx:498
AliAnalysisTaskPhiFlow::fV0Data
TProfile * fV0Data[18][2]
subevent resolution info for v2
Definition: AliAnalysisTaskPhiFlow.h:231
AliAnalysisTaskPhiFlow::fCentralityMin
Double_t fCentralityMin
QA profile of centralty vs multiplicity.
Definition: AliAnalysisTaskPhiFlow.h:207
Double_t
double Double_t
Definition: External.C:58
AliAnalysisTaskPhiFlow::CheckVertex
Bool_t CheckVertex(const T *event)
Definition: AliAnalysisTaskPhiFlow.cxx:382
AliFlowTrackSimple::SetEta
void SetEta(Double_t eta)
Definition: AliFlowTrackSimple.h:174
AliAnalysisTaskPhiFlow::fDCASecondaryWeak
TH2F * fDCASecondaryWeak
dca of primaries (mc) or kaons (data)
Definition: AliAnalysisTaskPhiFlow.h:228
AliAnalysisTaskPhiFlow::fPtKP
TH1F * fPtKP
QA histogram of p_t distribution of negative particles.
Definition: AliAnalysisTaskPhiFlow.h:203
TH2F
Definition: External.C:236
AliAnalysisTaskPhiFlow::fFlowEvent
AliFlowEvent * fFlowEvent
pid response object
Definition: AliAnalysisTaskPhiFlow.h:167
AliAnalysisTaskVnV0.h
AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing
void ReconstructionWithEventMixing(TObjArray *MixingCandidates) const
Definition: AliAnalysisTaskPhiFlow.cxx:972
AliAnalysisTaskPhiFlow::PhiPt
Double_t PhiPt(const T *track_1, const T *track_2) const
Definition: AliAnalysisTaskPhiFlow.cxx:681
AliFlowEventSimple::GetNumberOfRPs
Int_t GetNumberOfRPs() const
Definition: AliFlowEventSimple.h:59
AliFlowCommonConstants.h
AliFlowTrackCuts.h
AliAnalysisTaskPhiFlow::fPIDk
TH2F * fPIDk
QA histogram of TPC response of all charged particles.
Definition: AliAnalysisTaskPhiFlow.h:194
mass
Double_t mass
Definition: charmCutsOptimization.C:35
AliAnalysisTaskPhiFlow::fPtKN
TH1F * fPtKN
QA histogram of p_t distribution of positive kaons.
Definition: AliAnalysisTaskPhiFlow.h:204
AliAnalysisTaskPhiFlow::fDCAPrim
TH2F * fDCAPrim
qa plot of dca z
Definition: AliAnalysisTaskPhiFlow.h:227
AliFlowCandidateTrack::AddDaughter
void AddDaughter(Int_t value)
Definition: AliFlowCandidateTrack.h:27
AliAnalysisTaskPhiFlow::PlotMultiplcities
void PlotMultiplcities(const T *event) const
Definition: AliAnalysisTaskPhiFlow.cxx:373
AliAnalysisTaskPhiFlow::fDCAMaterial
TH2F * fDCAMaterial
dca of weak (mc)
Definition: AliAnalysisTaskPhiFlow.h:229
AliFlowBayesianPID::ComputeProb
void ComputeProb(const AliESDtrack *t, Float_t)
Definition: AliFlowBayesianPID.cxx:872
c
TCanvas * c
Definition: TestFitELoss.C:172
AliAnalysisTaskPhiFlow::fPtSpectra
TH1F * fPtSpectra[18]
like-sign kaon pairs
Definition: AliAnalysisTaskPhiFlow.h:200
AliFlowTrackSimple::SetMass
void SetMass(Double_t mass)
Definition: AliFlowTrackSimple.h:184
AliAnalysisTaskPhiFlow::fVZEROC
TH1F * fVZEROC
QA plot vzeroa multiplicity (all tracks in event)
Definition: AliAnalysisTaskPhiFlow.h:219
AliAnalysisTaskPhiFlow::InvariantMass
Double_t InvariantMass(const T *track1, const T *track2) const
Definition: AliAnalysisTaskPhiFlow.cxx:315
AliFlowTrackCuts::IsSelected
virtual Bool_t IsSelected(TObject *obj, Int_t id=-666)
Definition: AliFlowTrackCuts.cxx:734
AliFlowTrackCuts::SetParamType
void SetParamType(trackParameterType paramType)
Definition: AliFlowTrackCuts.h:227
AliFlowBayesianPID::SetDetResponse
void SetDetResponse(AliESDEvent *esd, Float_t centrality=-1.0, EStartTimeType_t flagStart=AliESDpid::kTOF_T0, Bool_t=kFALSE)
Definition: AliFlowBayesianPID.cxx:150
AliAnalysisTaskPhiFlow::fPIDTOF
TH2F * fPIDTOF
QA histo of TOF repsonse charged particles.
Definition: AliAnalysisTaskPhiFlow.h:196
AliFlowEventSimple::SetReferenceMultiplicity
void SetReferenceMultiplicity(Int_t m)
Definition: AliFlowEventSimple.h:54
AliAnalysisTaskVnV0::GetPsi2TPC
static Float_t GetPsi2TPC()
Definition: AliAnalysisTaskVnV0.h:61
AliFlowTrackSimple::Mass
Double_t Mass() const
Definition: AliFlowTrackSimple.h:161
AliFlowEventSimple.h
AliFlowTrackCuts::fV0
AliESDv0 * fV0
placeholder for the current kink
Definition: AliFlowTrackCuts.h:437
AliAnalysisTaskPhiFlow::GetWDist
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
Definition: AliAnalysisTaskPhiFlow.cxx:1151
AliFlowTrackSimple::InRPSelection
Bool_t InRPSelection() const
Definition: AliFlowTrackSimple.h:164
AliFlowTrackCuts::fBayesianResponse
AliFlowBayesianPID * fBayesianResponse
Definition: AliFlowTrackCuts.h:452
AliFlowEvent::GetTrack
AliFlowTrack * GetTrack(Int_t i)
Definition: AliFlowEvent.cxx:166
T
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
AliFlowBayesianPID::SetMC
void SetMC(Bool_t flag)
Definition: AliFlowBayesianPID.h:90
AliAnalysisTaskPhiFlow::fSkipEventSelection
Bool_t fSkipEventSelection
profiles for vzero vn(minv)
Definition: AliAnalysisTaskPhiFlow.h:232
AliFlowTrackSimple::SetForRPSelection
void SetForRPSelection(Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:188
AliAnalysisTaskPhiFlow::fPt
TH1F * fPt
QA plot of azimuthal distribution of tracks used for event plane estimation.
Definition: AliAnalysisTaskPhiFlow.h:216
AliFlowTrackCuts::SetPtRange
void SetPtRange(Float_t r1, Float_t r2)
Definition: AliFlowTrackCuts.h:135
AliFlowCandidateTrack::GetNDaughters
Int_t GetNDaughters(void) const
Definition: AliFlowCandidateTrack.h:26
AliFlowTrackSimple::Phi
Double_t Phi() const
Definition: AliFlowTrackSimple.h:155
AliAnalysisTaskPhiFlow::Exec
virtual void Exec(Option_t *)
Definition: AliAnalysisTaskPhiFlow.cxx:964
AliFlowEventSimple::SetNumberOfRPs
void SetNumberOfRPs(Int_t nr)
Definition: AliFlowEventSimple.h:60
AliFlowEvent::Fill
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
Definition: AliFlowEvent.cxx:858
Int_t
int Int_t
Definition: External.C:63
AliAnalysisTaskPhiFlow::PassesDCACut
Bool_t PassesDCACut(AliAODTrack *track) const
Definition: AliAnalysisTaskPhiFlow.cxx:512
AliAnalysisTaskPhiFlow::PrepareFlowEvent
void PrepareFlowEvent(Int_t iMulti)
Definition: AliAnalysisTaskPhiFlow.cxx:736
Float_t
float Float_t
Definition: External.C:68
AliFlowBayesianPID::GetCurrentMask
Bool_t GetCurrentMask(Int_t idet) const
Definition: AliFlowBayesianPID.h:106
AliFlowBayesianPID.h
AliAnalysisTaskPhiFlow::fVZEROA
TH1F * fVZEROA
QA plot of eta distribution of tracks used for event plane estimation.
Definition: AliAnalysisTaskPhiFlow.h:218
AliFlowCommonConstants::GetMaster
static AliFlowCommonConstants * GetMaster()
Definition: AliFlowCommonConstants.cxx:70
AliAnalysisTaskPhiFlow::fMultvsCentr
TH2F * fMultvsCentr
QA profile global and tpc multiplicity after outlier cut.
Definition: AliAnalysisTaskPhiFlow.h:206
AliAnalysisTaskVnV0::GetPsi2V0C
static Float_t GetPsi2V0C()
Definition: AliAnalysisTaskVnV0.h:60
AliAnalysisTaskPhiFlow::fEta
TH1F * fEta
QA plot of p_t sectrum of tracks used for event plane estimation.
Definition: AliAnalysisTaskPhiFlow.h:217
AliAnalysisTaskVnV0::GetPsi2V0A
static Float_t GetPsi2V0A()
Definition: AliAnalysisTaskVnV0.h:59
AliAnalysisTaskPhiFlow::fInvMNN
TH1F * fInvMNN[18]
unlike sign kaon pairs
Definition: AliAnalysisTaskPhiFlow.h:198
kProton
Definition: FastCentEstimators.C:23
AliAnalysisTaskPhiFlow::InitializeEventMixing
AliEventPoolManager * InitializeEventMixing()
Definition: AliAnalysisTaskPhiFlow.cxx:294
AliAnalysisTaskPhiFlow::PhiTrack
Bool_t PhiTrack(T *track) const
Definition: AliAnalysisTaskPhiFlow.cxx:718
AliFlowTrackSimple.h
AliFlowTrackCuts::fPIDResponse
AliPIDResponse * fPIDResponse
Definition: AliFlowTrackCuts.h:485
AliFlowTrackSimple::SetPhi
void SetPhi(Double_t phi)
Definition: AliFlowTrackSimple.h:178
AliAnalysisTaskPhiFlow::fMultCorAfterCuts
TH2F * fMultCorAfterCuts
QA histogram of p_t distribution of negative kaons.
Definition: AliAnalysisTaskPhiFlow.h:205
Short_t
short Short_t
Definition: External.C:23
AliFlowEventSimple::DefineDeadZone
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
Definition: AliFlowEventSimple.cxx:1139
AliFlowTrackSimple::SetPt
void SetPt(Double_t pt)
Definition: AliFlowTrackSimple.h:176
AliAnalysisTaskPhiFlow::fDCAZQA
TH1F * fDCAZQA
qa plot of dca xz
Definition: AliAnalysisTaskPhiFlow.h:226
AliFlowCandidateTrack::GetIDDaughter
Int_t GetIDDaughter(Int_t value) const
Definition: AliFlowCandidateTrack.h:28
AliAnalysisTaskPhiFlow::fCandidates
TObjArray * fCandidates
PID response object.
Definition: AliAnalysisTaskPhiFlow.h:169
AliAnalysisTaskPhiFlow::fPOICuts
AliFlowTrackCuts * fPOICuts
QA profile of centralty vs multiplicity.
Definition: AliAnalysisTaskPhiFlow.h:213
ClassImp
ClassImp(AliAnalysisTaskPhiFlow) ClassImp(AliPhiMesonHelperTrack) AliAnalysisTaskPhiFlow
Definition: AliAnalysisTaskPhiFlow.cxx:61
AliAnalysisTaskPhiFlow::fInvMNP
TH1F * fInvMNP[18]
QA histo of TOF response kaons.
Definition: AliAnalysisTaskPhiFlow.h:197
AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut
Bool_t CheckCandidateEtaPtCut(const T *track1, const T *track2) const
Definition: AliAnalysisTaskPhiFlow.cxx:349
AliFlowTrackCuts::SetEvent
void SetEvent(AliVEvent *event, AliMCEvent *mcEvent=NULL)
Definition: AliFlowTrackCuts.cxx:670
AliAnalysisTaskPhiFlow::InitializeBayesianPID
void InitializeBayesianPID(AliAODEvent *event)
Definition: AliAnalysisTaskPhiFlow.cxx:491
nEvents
Float_t nEvents[nProd]
Definition: DrawProductionComparison.C:36
AliAnalysisTaskPhiFlow::fPtN
TH1F * fPtN
QA histogram of p_t distribution of positive particles.
Definition: AliAnalysisTaskPhiFlow.h:202
AliAnalysisTaskPhiFlow::BookHistogram
TH1F * BookHistogram(const char *name)
Definition: AliAnalysisTaskPhiFlow.cxx:115
AliAnalysisTaskPhiFlow::fBayesianResponse
AliFlowBayesianPID * fBayesianResponse
flow events (one for each inv mass band)
Definition: AliAnalysisTaskPhiFlow.h:168
AliFlowTrackSimple::Pt
Double_t Pt() const
Definition: AliFlowTrackSimple.h:153
AliAnalysisTaskPhiFlow::GetDoubleCountingK
Bool_t GetDoubleCountingK(Double_t nSk, Short_t minNSigma) const
Definition: AliAnalysisTaskPhiFlow.cxx:1251
AliFlowBayesianPID::SetNewTrackParam
void SetNewTrackParam(Bool_t flag=kTRUE)
Definition: AliFlowBayesianPID.h:84
AliAnalysisTaskPhiFlow::InitPtSpectraHistograms
TH1F * InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
Definition: AliAnalysisTaskPhiFlow.cxx:146
AliFlowTrackSimple::SetForPOISelection
void SetForPOISelection(Bool_t b=kTRUE)
Definition: AliFlowTrackSimple.h:190
AliAnalysisTaskPhiFlow::fDCAXYQA
TH1F * fDCAXYQA
qa dca of all charged particles
Definition: AliAnalysisTaskPhiFlow.h:225
AliPhiMesonHelperTrack::Charge
virtual Int_t Charge() const
Definition: AliAnalysisTaskPhiFlow.h:39
AliFlowCandidateTrack.h
AliAnalysisTaskPhiFlow::fNOPIDTOF
TH2F * fNOPIDTOF
QA histogram of TPC response of kaons.
Definition: AliAnalysisTaskPhiFlow.h:195
Option_t
const char Option_t
Definition: External.C:48
AliAnalysisTaskPhiFlow::IsKaon
Bool_t IsKaon(AliAODTrack *track) const
Definition: AliAnalysisTaskPhiFlow.cxx:553
AliAnalysisTaskVnV0::IsPsiComputed
static Bool_t IsPsiComputed()
Definition: AliAnalysisTaskVnV0.h:58
AliAnalysisTaskPhiFlow::fDCAAll
TH2F * fDCAAll
QA plot TPC multiplicity (tracks used for event plane estimation)
Definition: AliAnalysisTaskPhiFlow.h:224
AliAnalysisTaskPhiFlow::fOutputList
TList * fOutputList
event pool manager
Definition: AliAnalysisTaskPhiFlow.h:189
Bool_t
bool Bool_t
Definition: External.C:53
AliAnalysisTaskPhiFlow::Terminate
virtual void Terminate(Option_t *)
Definition: AliAnalysisTaskPhiFlow.cxx:1063
AliAnalysisTaskPhiFlow::fSubEventDPhiv2
TProfile * fSubEventDPhiv2
dca material (mc) all (data)
Definition: AliAnalysisTaskPhiFlow.h:230
AliAnalysisTaskPhiFlow::BookPtHistogram
TH1F * BookPtHistogram(const char *name)
Definition: AliAnalysisTaskPhiFlow.cxx:156
AliAnalysisTaskPhiFlow.h
AliAnalysisTaskPhiFlow::fTPCM
TH1F * fTPCM
QA plot vzeroc multiplicity (all tracks in event)
Definition: AliAnalysisTaskPhiFlow.h:220
AliAnalysisTaskPhiFlow::FindMinNSigma
Short_t FindMinNSigma(Double_t nSpi, Double_t nSk, Double_t nSp) const
Definition: AliAnalysisTaskPhiFlow.cxx:1230
AliAnalysisTaskPhiFlow::BookPIDHistogram
TH2F * BookPIDHistogram(const char *name, Bool_t TPC)
Definition: AliAnalysisTaskPhiFlow.cxx:128
AliAnalysisTaskPhiFlow::UserExec
virtual void UserExec(Option_t *option)
Definition: AliAnalysisTaskPhiFlow.cxx:799
AliFlowEvent::InsertTrack
void InsertTrack(AliFlowTrack *)
Definition: AliFlowEvent.cxx:1017
AliFlowEvent::ClearFast
virtual void ClearFast()
Definition: AliFlowEvent.cxx:2406
AliAnalysisTaskPhiFlow::fPoolManager
AliEventPoolManager * fPoolManager
AOD oject.
Definition: AliAnalysisTaskPhiFlow.h:188
AliAnalysisTaskPhiFlow::fInvMPP
TH1F * fInvMPP[18]
like-sign kaon pairs
Definition: AliAnalysisTaskPhiFlow.h:199
AliAnalysisTaskPhiFlow::PtSelector
void PtSelector(Int_t _track_type, const T *track_1, const T *track_2) const
Definition: AliAnalysisTaskPhiFlow.cxx:690
AliAnalysisTaskPhiFlow::fNOPID
TH2F * fNOPID
QA histogram of events that do not pass centrality cut.
Definition: AliAnalysisTaskPhiFlow.h:193
AliAnalysisTaskPhiFlow::MakeTrack
void MakeTrack(Double_t, Double_t, Double_t, Double_t, Int_t, Int_t[], Double_t p=0., Double_t pz=0.) const
Definition: AliAnalysisTaskPhiFlow.cxx:1069
AliFlowTrackCuts::SetEtaRange
void SetEtaRange(Float_t r1, Float_t r2)
Definition: AliFlowTrackCuts.h:134
AliFlowEventSimple::NumberOfTracks
Int_t NumberOfTracks() const
Definition: AliFlowEventSimple.h:52