AliPhysics  c923f52 (c923f52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskCascadeTester.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 // AliAnalysisTaskCascadeTester:
17 // author: Redmer Alexander Bertens (rbertens@cern.ch)
18 
19 #include "TChain.h"
20 #include "TH1.h"
21 #include "TH1F.h"
22 #include "TH2F.h"
23 #include "TMath.h"
24 #include "TObjArray.h"
25 #include "AliAnalysisTaskSE.h"
26 #include "AliAnalysisManager.h"
27 #include "AliAODEvent.h"
28 #include "AliAODInputHandler.h"
29 #include "AliCentrality.h"
31 #include "AliPIDCombined.h"
32 #include "AliMCEvent.h"
33 #include "TProfile.h"
34 #include "AliFlowCandidateTrack.h"
35 #include "AliFlowTrackCuts.h"
36 #include "AliFlowEventSimple.h"
37 #include "AliFlowTrackSimple.h"
38 #include "AliFlowCommonConstants.h"
39 #include "AliFlowEvent.h"
40 #include "TVector3.h"
41 #include "AliAODVZERO.h"
42 #include "AliPIDResponse.h"
43 #include "AliTPCPIDResponse.h"
44 #include "AliAODMCParticle.h"
45 #include "AliAnalysisTaskVnV0.h"
46 #include "AliEventPoolManager.h"
47 #include "AliMultSelection.h"
48 #include "TMatrixDSym.h"
49 #include "AliVVertex.h"
50 #include "AliAODcascade.h"
51 
52 class AliFlowTrackCuts;
53 
54 using std::cout;
55 using std::endl;
56 
58 
60  fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0),fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0), fPileUp(kTRUE), fOmega(kTRUE), fXi(kFALSE)
61 {
62  // Default constructor
63  for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 1000.;
64  for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
65  for(Int_t i(0); i < 20; i++) {
66  fVertexMixingBins[i] = 0;
67  fCentralityMixingBins[i] = 0;
68  }
69  fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
70  for(Int_t i(0); i < 18; i++) {
71  for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
72  fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
73  }
74 }
75 //_____________________________________________________________________________
77  fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fCandidateYCut(kFALSE), fCandidateMinY(-.5), fCandidateMaxY(.5), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0), fPileUp(kTRUE), fOmega(kTRUE), fXi(kFALSE)
78 {
79  // Constructor
80  for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 1000.;
81  for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
82  for(Int_t i(0); i < 20; i++) {
83  fVertexMixingBins[i] = 0;
84  fCentralityMixingBins[i] = 0;
85  }
86  fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
87  for(Int_t i(0); i < 18; i++) {
88  for(Int_t j(0); j < 2; j++) fV0Data[i][j] = 0;
89  fInvMNP[i] = 0; fInvMNN[i] = 0; fInvMPP[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.;
90  }
91  DefineInput(0, TChain::Class());
92  DefineOutput(1, TList::Class());
93  DefineOutput(2, AliFlowEventSimple::Class());
94  if(fDebug > 0) cout << " === Created instance of AliAnaysisTaskCascadeTester === " << endl;
95 }
96 //_____________________________________________________________________________
98 {
99  // Destructor
100  if (fNullCuts) delete fNullCuts;
101  if (fOutputList) delete fOutputList;
102  if (fCandidates) delete fCandidates;
103  if (fFlowEvent) delete fFlowEvent;
104  if (fEventMixing) delete fPoolManager;
105  if (fPIDCombined) delete fPIDCombined;
106  if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskCascadeTester === " << endl;
107 }
108 //_____________________________________________________________________________
110 {
111  // Return a pointer to a TH1 with predefined binning
112  if(fDebug > 0) cout << " *** BookHistogram() *** " << name << endl;
113  TH1F *hist = 0x0;
114  if(fOmega) hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, 1.64, 1.71);
115  else if(fXi) hist = new TH1F(name, Form("M_{INV} (%s)", name), 60, 1.29, 1.36);
116  hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
117  hist->GetYaxis()->SetTitle("No. of pairs");
118  hist->SetMarkerStyle(kFullCircle);
119  hist->Sumw2();
120  fOutputList->Add(hist);
121  return hist;
122 }
123 //_____________________________________________________________________________
125 {
126  // Return a pointer to a TH2 with predefined binning
127  if(fDebug > 0) cout << " *** BookPIDHisotgram() *** " << endl;
128  TH2F *hist = 0x0;
129  if(TPC) {
130  hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
131  hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
132  }
133  if(!TPC) {
134  hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
135  hist->GetYaxis()->SetTitle("#beta");
136  }
137  hist->GetXaxis()->SetTitle("P (GeV / c)");
138  fOutputList->Add(hist);
139  return hist;
140 }
141 //_____________________________________________________________________________
143 {
144  // intialize p_t histograms for each p_t bin
145  if(fDebug > 0) cout << " *** InitPtSpectraHistograms() *** " << endl;
146  TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 60, nmin, nmax);
147  hist->GetXaxis()->SetTitle("p_{T} GeV / c");
148  fOutputList->Add(hist);
149  return hist;
150 }
151 //_____________________________________________________________________________
153 {
154  // Return a pointer to a p_T spectrum histogram
155  if(fDebug > 0) cout << " *** BookPtHistogram() *** " << endl;
156  TH1F* ratio = new TH1F(name, name, 100, 0, 7);
157  ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
158  ratio->GetYaxis()->SetTitle("No. of events");
159  ratio->Sumw2();
160  fOutputList->Add(ratio);
161  return ratio;
162 }
163 //_____________________________________________________________________________
165 {
166  // Add Phi Identification Output Objects
167  if(fDebug > 0) cout << " ** AddPhiIdentificationOutputObjects() *** " << endl;
168  if(fQA) {
169  fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
170  fEventStats->GetXaxis()->SetTitle("No. of events");
171  fEventStats->GetYaxis()->SetTitle("Statistics");
172  fOutputList->Add(fEventStats);
173  }
174  fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
176  if(fQA) {
177  fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
179  fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
180  fPIDk = BookPIDHistogram("TPC signal, kaons", kTRUE);
181  fNOPIDTOF = BookPIDHistogram("TOF signal, all particles", kFALSE);
182  fPIDTOF = BookPIDHistogram("TOF signal, kaons", kFALSE);
183  }
184  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
185  fInvMNP[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
186  fInvMNN[ptbin] = BookHistogram(Form("NN, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
187  fInvMPP[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
188  }
189  if(fQA) {
190  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
191  fPtP = BookPtHistogram("i^{+}");
192  fPtN = BookPtHistogram("i^{-}");
193  fPtKP = BookPtHistogram("K^{+}");
194  fPtKN = BookPtHistogram("K^{-}");
195  fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
196  fOutputList->Add(fPhi);
197  fPt = new TH1F("fPt", "p_{T}", 100, 0, 5.5);
198  fOutputList->Add(fPt);
199  fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
200  fOutputList->Add(fEta);
201  fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
202  fOutputList->Add(fVZEROA);
203  fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
204  fOutputList->Add(fVZEROC);
205  fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
206  fOutputList->Add(fTPCM);
207  fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
208  fOutputList->Add(fDCAXYQA);
209  fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
210  fOutputList->Add(fDCAZQA);
212  fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
214  fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
216  }
217  }
218  if(fIsMC || fQA) {
219  fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
220  fOutputList->Add(fDCAAll);
221  fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
222  fOutputList->Add(fDCAPrim);
223  fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
225  fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
227  }
228  if(fV0) {
229  fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 5, 0, 5);
230  fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<#Psi_{a} - #Psi_{b}>");
231  fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<#Psi_{a} - #Psi_{c}>");
232  fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<#Psi_{b} - #Psi_{c}>");
233  fSubEventDPhiv2->GetXaxis()->SetBinLabel(4, "#sqrt{#frac{<#Psi_{a} - #Psi_{b}><#Psi_{a} - #Psi_{c}>}{<#Psi_{b} - #Psi_{c}>}}");
234  fSubEventDPhiv2->GetXaxis()->SetBinLabel(5, "#sqrt{#frac{<#Psi_{a} - #Psi_{c}><#Psi_{b} - #Psi_{c}>}{<#Psi_{a} - #Psi_{b}>}}");
236  const char* V0[] = {"V0A", "V0C"};
237  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
238  for(Int_t iV0(0); iV0 < 2; iV0++) {
239  fV0Data[ptbin][iV0] = new TProfile(Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), Form("%s v2 %4.2f < p_{T} < %4.2f GeV", V0[iV0], fPtBins[ptbin], fPtBins[ptbin+1]), fMassBins, fMinMass, fMaxMass);
240  fOutputList->Add(fV0Data[ptbin][iV0]);
241  }
242  }
243 }
244 //_____________________________________________________________________________
246 {
247  // Create user defined output objects
248  if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
249  fNullCuts = new AliFlowTrackCuts("null_cuts");
250  // combined pid
251  fPIDCombined = new AliPIDCombined;
252  fPIDCombined->SetDefaultTPCPriors();
253  fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
254 
255  Double_t t(0);
256  for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
257  if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
258  if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
260  cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
261  cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
262  cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
263  cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
264  cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
265  cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
266  AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
267  if (man) {
268  AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
269  if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
270  }
271  // Create all output objects and store them to a list
272  fOutputList = new TList();
273  fOutputList->SetOwner(kTRUE);
274  // Create and post the Phi identification output objects
276  PostData(1, fOutputList);
277  // create candidate array
278  fCandidates = new TObjArray(1000);
279  fCandidates->SetOwner(kTRUE);
280  // create and post flowevent
281  fFlowEvent = new AliFlowEvent(10000);
282  PostData(2, fFlowEvent);
284 }
285 //_____________________________________________________________________________
287 {
288  // initialize event mixing
289  if(fDebug > 0) cout << " *** InitializeEventMixing() *** " << endl;
290  Int_t _c(0), _v(0);
291  for(Int_t i(0); i < 19; i++) {
292  if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
293  else _c = 19;
294  }
295  for(Int_t i(0); i < 19; i++) {
296  if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
297  else _v = 19;
298  }
299  if(fDebug > 0 ) cout << Form(" --> found %d centrality bins for mixing, %d vertex bins for mixing", _c, _v) << endl;
300  Double_t centralityBins[_c];
301  Double_t vertexBins[_v];
302  for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
303  for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
304  return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
305 }
306 //_____________________________________________________________________________
307 template <typename T> Double_t AliAnalysisTaskCascadeTester::InvariantMass(const T* track1, const T* track2) const
308 {
309  // Return the invariant mass of two tracks, assuming both tracks are kaons
310  if(fDebug > 1) cout << " *** InvariantMass() *** " << endl;
311  if ((!track2) || (!track1)) return 0.;
312  Double_t masss = TMath::Power(4.93676999999999977e-01, 2);
313  Double_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
314  Double_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
315  Double_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
316  Double_t e1 = TMath::Sqrt(track1->P() * track1->P() + masss);
317  Double_t e2 = TMath::Sqrt(track2->P() * track2->P() + masss);
318  Double_t es = TMath::Power((e1 + e2), 2);
319  if ((es - (pxs + pys + pzs)) < 0) return 0.;
320  return TMath::Sqrt((es - (pxs + pys + pzs)));
321 }
322 //_____________________________________________________________________________
323 /*
324 template <typename T> Double_t AliAnalysisTaskCascadeTester::DeltaDipAngle(const T* track1, const T* track2) const
325 {
326  // Calculate the delta dip angle between two particles
327  if(fDebug > 1) cout << " *** DeltaDipAngle() *** " << endl;
328  if (track1->P()*track2->P() == 0) return 999;
329  return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
330 }
331 //_____________________________________________________________________________
332 template <typename T> Bool_t AliAnalysisTaskCascadeTester::CheckDeltaDipAngle(const T* track1, const T* track2) const
333 {
334  // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
335  if(fDebug > 1) cout << " *** CheckDeltaDipAngle() *** " << endl;
336  if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
337  return kTRUE;
338 }
339 */
340 //_____________________________________________________________________________
341 template <typename T> Bool_t AliAnalysisTaskCascadeTester::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
342 {
343  // Check if pair passes eta and pt cut
344  if(fDebug > 1) cout << " *** CheckCandidateEtaPtCut() *** " << endl;
345  if (fCandidateMinPt > PhiPt(track1, track2) || fCandidateMaxPt < PhiPt(track1, track2)) return kFALSE;
346  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
347  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
348  TVector3 c = a + b;
349  if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
350  return kTRUE;
351 }
352 //_____________________________________________________________________________
353 template <typename T> Bool_t AliAnalysisTaskCascadeTester::EventCut(T* event)
354 {
355  // Impose event cuts
356  if(fDebug > 0) cout << " *** EventCut() *** " << endl;
357  if (!event) return kFALSE;
358 // if (fSkipEventSelection) return kTRUE;
359  if (!CheckVertex(event)) return kFALSE;
360  if (!CheckCentrality(event)) return kFALSE;
361  if(fQA) PlotMultiplcities(event);
362  return kTRUE;
363 }
364 //_____________________________________________________________________________
365 template <typename T> void AliAnalysisTaskCascadeTester::PlotMultiplcities(const T* event) const
366 {
367  // QA multiplicity plots
368  if(fDebug > 1) cout << " *** PlotMultiplcities() *** " << endl;
369  fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
370  fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
371  fTPCM->Fill(event->GetNumberOfTracks());
372 }
373 //_____________________________________________________________________________
374 template <typename T> Bool_t AliAnalysisTaskCascadeTester::CheckVertex(const T* event)
375 {
376  // Check if event vertex is within given range
377  if(fDebug > 0) cout << " *** CheckVertex() *** " << endl;
378  if (!event->GetPrimaryVertex()) return 0x0;
379  fVertex = event->GetPrimaryVertex()->GetZ();
380  if (TMath::Abs(fVertex) > fVertexRange) return 0x0;
381  return kTRUE;
382 }
383 //_____________________________________________________________________________
385 {
386  // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
387  if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
388  if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
389 
390  // check if the AliMultSelection object is present. If so, we should invoke the
391  // new centrality framework
392 
393  AliMultSelection *multSelection = 0x0;
394  multSelection = static_cast<AliMultSelection*>(event->FindListObject("MultSelection"));
395  if(multSelection) {
396  fCentrality = multSelection->GetMultiplicityPercentile("V0M");
399  return kTRUE;
400  } else {
402  return kFALSE;
403  }
404  }
405 
406  else fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
407  Double_t cenB(-999);
408  // if a second centrality estimator is requited, set it
409  (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
410  if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
411  if(fQA) fCentralityNoPass->Fill(fCentrality) ;
412  return kFALSE;
413  }
414  const Int_t nGoodTracks = event->GetNumberOfTracks();
415  if(fCentralityCut2010) { // cut on outliers
416  Float_t multTPC(0.); // tpc mult estimate
417  Float_t multGlob(0.); // global multiplicity
418  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
419  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
420  if(!trackAOD) AliFatal("Not a standard AOD");
421  if (!trackAOD) continue;
422  if (!(trackAOD->TestFilterBit(1))) continue;
423  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;
424  multTPC++;
425  }
426  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global 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(16))) 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.1)) continue;
432  Double_t b[2] = {-99., -99.};
433  Double_t bCov[3] = {-99., -99., -99.};
434  AliAODTrack copy(*trackAOD);
435  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
436  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
437  multGlob++;
438  } //track loop
439  // printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
440  if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
441  if(fQA) {
442  fMultCorAfterCuts->Fill(multGlob, multTPC);
443  fMultvsCentr->Fill(fCentrality, multTPC);
444  }
445  }
446 
447  if(fCentralityCut2011) { // cut on outliers
448  Float_t multTPC(0.); // tpc mult estimate
449  Float_t multGlob(0.); // global multiplicity
450  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
451  AliAODTrack* trackAOD = dynamic_cast<AliAODTrack*>(event->GetTrack(iTracks));
452  if(!trackAOD) AliFatal("Not a standard AOD");
453  if (!trackAOD) continue;
454  if (!(trackAOD->TestFilterBit(1))) continue;
455  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;
456  multTPC++;
457  }
458  for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global 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(16))) 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.1)) continue;
464  Double_t b[2] = {-99., -99.};
465  Double_t bCov[3] = {-99., -99., -99.};
466  AliAODTrack copy(*trackAOD);
467  if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
468  if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
469  multGlob++;
470  } //track loop
471  //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
472  if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
473  if(fQA) {
474  fMultCorAfterCuts->Fill(multGlob, multTPC);
475  fMultvsCentr->Fill(fCentrality, multTPC);
476  }
477  }
478 
480  return kTRUE;
481 }
482 //_____________________________________________________________________________
483 template <typename T> Double_t AliAnalysisTaskCascadeTester::PhiPt(const T* track1, const T* track2) const
484 {
485  // return p_t of track pair
486  TVector3 a(track1->Px(), track1->Py(), track1->Pz());
487  TVector3 b(track2->Px(), track2->Py(), track2->Pz());
488  TVector3 c = a + b;
489  return c.Pt();
490 }
491 //_____________________________________________________________________________
492 template <typename T> void AliAnalysisTaskCascadeTester::PtSelector(Int_t tracktype, const T* track1, const T* track2) const
493 {
494  // plot m_inv spectra of like- and unlike-sign kaon pairs for each pt bin
495  Double_t pt = PhiPt(track1, track2);
496  if (tracktype == 0) {
497  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
498  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
499  fInvMNP[ptbin]->Fill(InvariantMass(track1, track2));
500  if(fQA) fPtSpectra[ptbin]->Fill(pt);
501  }
502  }
503  }
504  if (tracktype == 1) {
505  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
506  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
507  fInvMPP[ptbin]->Fill(InvariantMass(track1, track2));
508  }
509  }
510  }
511  if (tracktype == 2) {
512  for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
513  if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
514  fInvMNN[ptbin]->Fill(InvariantMass(track1, track2));
515  }
516  }
517  }
518 }
519 //_____________________________________________________________________________
520 template <typename T> Bool_t AliAnalysisTaskCascadeTester::PhiTrack(T* track) const
521 {
522  // check if track meets quality cuts
523  if(!track) return kFALSE;
524  return fPOICuts->IsSelected(track);
525 }
526 //_____________________________________________________________________________
527 template <typename T> void AliAnalysisTaskCascadeTester::SetNullCuts(T* event)
528 {
529  // Set null cuts
530  if (fDebug > 0) cout << " *** SetNullCuts() *** " << fCutsRP << endl;
531  fCutsRP->SetEvent(event, MCEvent());
533  fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
534  fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
535  fNullCuts->SetEvent(event, MCEvent());
536 }
537 //_____________________________________________________________________________
539 {
540  // Prepare flow events
541  if (fDebug > 0 ) cout << " *** PrepareFlowEvent() *** " << endl;
545  fFlowEvent->DefineDeadZone(0, 0, 0, 0);
546 }
547 //_____________________________________________________________________________
549 {
550  // UserExec: called for each event. Commented where necessary
551  if(fDebug > 0 ) cout << " *** UserExec() *** " << endl;
552  TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
553  if(fEventMixing) {
554  MixingCandidates = new TObjArray();
555  MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
556  }
557  if (!fPIDResponse) {
558  if(fDebug > 0 ) cout << " Could not get PID response " << endl;
559  return;
560  }
561  fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
562  if (fAOD) {
563  if (!EventCut(fAOD)) return; // check for event cuts
564 
565 
566  // add extra pileup cuts for high intensity runs
567  // courtesy of alex dobrin
568  if (fPileUp){
569 
570  if (plpMV(fAOD))
571  return;
572 
573  Short_t isPileup = fAOD->IsPileupFromSPD(3);
574  if (isPileup != 0)
575  return;
576 
577  if (((AliAODHeader*)fAOD->GetHeader())->GetRefMultiplicityComb08() < 0)
578  return;
579 
580  /*
581  Int_t bc2 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt2ClosestInteractionMap();
582  if (bc2 != 0)
583  return;
584 
585  Int_t bc1 = ((AliAODHeader*)fAOD->GetHeader())->GetIRInt1ClosestInteractionMap();
586  if (bc1 != 0)
587  return;
588  */
589 
590  }
591 
592  //new function for 2015 to remove incomplete events
593  if (fAOD->IsIncompleteDAQ())
594  return;
595 
596 
597  SetNullCuts(fAOD);
598  PrepareFlowEvent(fAOD->GetNumberOfTracks());
599  fCandidates->SetLast(-1);
600  if(fIsMC) IsMC(); // launch mc stuff
601  if(fQA) fEventStats->Fill(0);
602 
603 
604  const Int_t ncasc = fAOD->GetNumberOfCascades();
605  Double_t lPrimaryVtxPosition[3];
606  fAOD->GetPrimaryVertex()->GetXYZ(lPrimaryVtxPosition);
607 
608  //cascades vn
609  for (Int_t jCasc = 0; jCasc < ncasc; jCasc++) {
610 
611  AliAODcascade* casc = (AliAODcascade*) fAOD->GetCascade(jCasc);
612 
613  if (!casc){
614  delete casc;
615  continue;
616  }
617 
618  Double_t lDcaXiDaughters = casc->DcaXiDaughters();
619  Double_t lDcaBachToPrimVertexXi = casc->DcaBachToPrimVertex();
620  Double_t lXiCosineOfPointingAngle = casc->CosPointingAngleXi(lPrimaryVtxPosition[0], lPrimaryVtxPosition[1], lPrimaryVtxPosition[2]);
621 
622  Double_t lPosXi[3] = { -1000.0, -1000.0, -1000.0 };
623  lPosXi[0] = casc->DecayVertexXiX();
624  lPosXi[1] = casc->DecayVertexXiY();
625  lPosXi[2] = casc->DecayVertexXiZ();
626 
627  Double_t lDcaV0DaughtersXi = casc->DcaV0Daughters();
628  Double_t lV0toXiCosineOfPointingAngle = casc->CosPointingAngle(lPosXi);
629 
630  Double_t lPosV0Xi[3] = { -1000.0, -1000.0, -1000.0 };
631  lPosV0Xi[0] = casc->DecayVertexV0X();
632  lPosV0Xi[1] = casc->DecayVertexV0Y();
633  lPosV0Xi[2] = casc->DecayVertexV0Z();
634 
635  Double_t lDcaV0ToPrimVertexXi = casc->DcaV0ToPrimVertex();
636  Double_t lDcaPosToPrimVertexXi = casc->DcaPosToPrimVertex();
637  Double_t lDcaNegToPrimVertexXi = casc->DcaNegToPrimVertex();
638 
639  Double_t lInvMassLambdaAsCascDghter = casc->MassLambda();
640 
641  Double_t lXiRadius = TMath::Sqrt(lPosXi[0]*lPosXi[0] + lPosXi[1]*lPosXi[1]);
642  Double_t lV0RadiusXi = TMath::Sqrt(lPosV0Xi[0]*lPosV0Xi[0] + lPosV0Xi[1]*lPosV0Xi[1]);
643 
644 
645  if (lDcaXiDaughters > 0.3)
646  continue; // in AliCascadeVertexer: max allowed DCA between the V0 and the bachelor
647 
648  if (lXiCosineOfPointingAngle < 0.999 )
649  continue; // in AliCascadeVertexer: min allowed cosine of the cascade pointing angle
650 
651  if (lDcaV0ToPrimVertexXi < 0.05)
652  continue; // in AliCascadeVertexer: min allowed V0 impact parameter
653 
654  if (lDcaBachToPrimVertexXi < 0.03)
655  continue; // in AliCascadeVertexer: min allowed bachelor's impact parameter
656 
657  if (TMath::Abs(lInvMassLambdaAsCascDghter-1.11568) > 0.008 )
658  continue; // in AliCascadeVertexer: "window" around the Lambda mass
659 
660  if(lXiRadius < 0.9)
661  continue; // in AliCascadeVertexer: min radius of the fiducial volume
662 
663  if(lXiRadius > 100.)
664  continue; // in AliCascadeVertexer: max radius of the fiducial volume
665 
666 
667  if (lDcaV0DaughtersXi > 1.)
668  continue; // in AliV0vertexer: max allowed DCA between the daughter tracks
669 
670  if (lV0toXiCosineOfPointingAngle < 0.998)
671  continue; // in AliV0vertexer: min allowed cosine of V0's pointing angle
672 
673  if (lDcaPosToPrimVertexXi < 0.1)
674  continue; // in AliV0vertexer: min allowed impact parameter for the 1st daughter
675 
676  if (lDcaNegToPrimVertexXi < 0.1)
677  continue; // in AliV0vertexer: min allowed impact parameter for the 2nd daughter
678 
679  if(lV0RadiusXi < 0.9)
680  continue; // in AliV0vertexer: min radius of the fiducial volume
681 
682  if(lV0RadiusXi > 100.)
683  continue; // in AliV0vertexer: max radius of the fiducial volume
684 
685 
686  if ((casc->Pt() < 0.2) || (casc->Pt() >= 8.) || (TMath::Abs(casc->Eta()) >= 0.8))
687  continue;
688 
689 
690 
691  AliAODTrack* pTrackXi = (AliAODTrack*)casc->GetDaughter(0);
692  AliAODTrack* nTrackXi = (AliAODTrack*)casc->GetDaughter(1);
693  AliAODTrack* bachTrackXi = (AliAODTrack*)casc->GetDecayVertexXi()->GetDaughter(0);
694 
695  if (!pTrackXi || !nTrackXi || !bachTrackXi)
696  continue;
697 
698 
699  UInt_t lIdxPosXi = (UInt_t) TMath::Abs(pTrackXi->GetID());
700  UInt_t lIdxNegXi = (UInt_t) TMath::Abs(nTrackXi->GetID());
701  UInt_t lBachIdx = (UInt_t) TMath::Abs(bachTrackXi->GetID());
702  if(lBachIdx == lIdxNegXi)
703  continue;
704 
705  if(lBachIdx == lIdxPosXi)
706  continue;
707 
708 
709  if (!(pTrackXi->IsOn(AliAODTrack::kTPCrefit)))
710  continue;
711 
712  if (!(nTrackXi->IsOn(AliAODTrack::kTPCrefit)))
713  continue;
714 
715  if (!(bachTrackXi->IsOn(AliAODTrack::kTPCrefit)))
716  continue;
717 
718 
719  if (pTrackXi->GetTPCNcls() < 70 || nTrackXi->GetTPCNcls() < 70 || bachTrackXi->GetTPCNcls() < 70)
720  continue;
721 
722  if (TMath::Abs(pTrackXi->Eta()) > 0.8 || (TMath::Abs(nTrackXi->Eta()) > 0.8 || TMath::Abs(bachTrackXi->Eta()) > 0.8))
723  continue;
724 
725  if (pTrackXi->Pt() < 0.15 || nTrackXi->Pt() < 0.15 || bachTrackXi->Pt() < 0.15)
726  continue;
727 
728 
729  //PID cuts
730  if (bachTrackXi->GetTPCsignalN() < 70 || nTrackXi->GetTPCsignalN() < 70 || pTrackXi->GetTPCsignalN() < 70)
731  continue;
732 
733  // Bachelor
734  Double_t nSigmaBachK = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bachTrackXi, AliPID::kKaon));
735  Double_t nSigmaBachPi = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(bachTrackXi, AliPID::kPion));
736 
737  // Negative V0 daughter
738  Double_t nSigmaPiNeg = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrackXi, AliPID::kPion));
739  Double_t nSigmaPNeg = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(nTrackXi,AliPID::kProton));
740 
741  // Positive V0 daughter
742  Double_t nSigmaPiPos = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrackXi,AliPID::kPion));
743  Double_t nSigmaPPos = TMath::Abs(fPIDResponse->NumberOfSigmasTPC(pTrackXi,AliPID::kProton));
744 
745 
746  Double_t lInvMassOmega = casc->MassOmega();
747  Double_t lInvMassXi = casc->MassXi();
748 
749  if (lInvMassOmega > 1.64 && lInvMassOmega < 1.71 && nSigmaBachK < 3.){
750  if ((nSigmaPiNeg < 3. && nSigmaPPos < 3.) || (nSigmaPNeg < 3. && nSigmaPiPos < 3.)){
751  if(fOmega) MakeTrack(lInvMassOmega, casc->Pt(), casc->Phi(), casc->Eta());
752  }
753 
754  }
755 
756  if (lInvMassXi > 1.29 && lInvMassXi < 1.36 && nSigmaBachPi < 3.){
757 
758  if ((nSigmaPiNeg < 3. && nSigmaPPos < 3.) || (nSigmaPNeg < 3. && nSigmaPiPos < 3.)){
759 
760  if(fXi) MakeTrack(lInvMassXi,casc->Pt(), casc->Phi(), casc->Eta());
761  //fill Xi invariant mass and v2 vs invariant mass hists
762  }
763 
764  }
765  }
766 
767 
768 
769  if (fDebug > 0) printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates into flow events
770  for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
771  AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
772  if (!cand) continue;
773  if (fDebug > 1) printf(" --> Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
774  for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
775  if (fDebug>1) printf(" *** Daughter %d with fID %d ***", iDau, cand->GetIDDaughter(iDau));
776  for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
777  AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
778  if (!iRP) continue;
779  if (!iRP->InRPSelection()) continue;
780  if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
781  if (fDebug > 1) printf(" was in RP set");
782  iRP->SetForRPSelection(kFALSE);
784  }
785  }
786  if (fDebug > 1) printf("\n");
787  }
788  cand->SetForPOISelection(kTRUE);
789  fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
790  }
791  if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
792 
793  /* if(!fEventMixing) { // combinatorial background
794  for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
795  for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
796 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
797  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
798  PtSelector(1, up[pTracks], up[nTracks]);
799  }
800  }
801  for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
802  for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
803 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
804  if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
805  PtSelector(2, un[nTracks], un[pTracks]);
806  }
807  }
808  }
809  if(fEventMixing) ReconstructionWithEventMixing(MixingCandidates);
810  */
811  PostData(1, fOutputList);
812  PostData(2, fFlowEvent);
813  }
814 }
815 //_____________________________________________________________________________
817 {
818  // skip the event selection for SE task (e.g. for MC productions)
820  // exec of task se will do event selection and call UserExec
821  else AliAnalysisTaskSE::Exec(c);
822 }
823 //_____________________________________________________________________________
825 {
826  // perform phi reconstruction with event mixing
827  if(fDebug > 0) cout << " *** ReconstructionWithEventMixing() *** " << endl;
828  AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
829  if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, fVertex));
830  if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
831  Int_t nEvents = pool->GetCurrentNEvents();
832  if(fDebug > 0) cout << " --> " << nEvents << " events in mixing buffer ... " << endl;
833  for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
834  TObjArray* mixed_candidates = pool->GetEvent(iEvent);
835  if(!mixed_candidates) continue; // this should NEVER happen
836  Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
837  Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
838  if(fDebug > 0) cout << Form(" - mixing %d tracks with %d buffered tracks from event %d ... ", candidates, bufferTracks, iEvent) << endl;
839  AliPhiMesonHelperTrack* buffer_un[bufferTracks]; // set values for buffered candidates
840  AliPhiMesonHelperTrack* buffer_up[bufferTracks];
841  Int_t buffer_unp(0);
842  Int_t buffer_unn(0);
843  AliPhiMesonHelperTrack* mix_un[candidates];// set values for mixing candidates
844  AliPhiMesonHelperTrack* mix_up[candidates];
845  Int_t mix_unp(0);
846  Int_t mix_unn(0);
847  for (Int_t iTracks = 0; iTracks < candidates; iTracks++) { // distinguish between k+ and k- for mixing candidates
848  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)MixingCandidates->At(iTracks);
849  if(!track) continue;
850  if (track->Charge() > 0) {
851  mix_up[mix_unp] = track;
852  mix_unp++;
853  }
854  else if (track->Charge() < 0 ) {
855  mix_un[mix_unn] = track;
856  mix_unn++;
857  }
858  }
859  for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) { // distinguish between k+ and k- for buffered candidates
860  AliPhiMesonHelperTrack* track = (AliPhiMesonHelperTrack*)mixed_candidates->At(iTracks);
861  if(!track) continue;
862  if (track->Charge() > 0) {
863  buffer_up[buffer_unp] = track;
864  buffer_unp++;
865  }
866  else if (track->Charge() < 0 ) {
867  buffer_un[buffer_unn] = track;
868  buffer_unn++;
869  }
870  }
871  for (Int_t pMix = 0; pMix < mix_unp; pMix++) { // mix k+ (current event) k+ (buffer)
872  if(fDebug > 1 ) cout << Form("mixing current k+ track %d with", pMix);
873  if(!fTypeMixing) { // mix with like-sign kaons
874  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
875  if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
876 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
877  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
878  PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
879  }
880  }
881  else if(fTypeMixing) { // mix with unlike-sign kaons
882  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
883  if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
884 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
885  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
886  PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
887  }
888  }
889  }
890  for (Int_t nMix = 0; nMix < mix_unn; nMix++) { // mix k- (current event) k- (buffer)
891  if(fDebug > 1 ) cout << Form("mixing current k- track %d with", nMix);
892  if(!fTypeMixing) { // mix with like-sign kaons
893  for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
894  if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
895 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
896  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
897  PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
898  }
899  }
900  else if(fTypeMixing) { // mix with unlike-sign kaons
901  for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
902  if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
903 // if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
904  if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
905  PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
906  }
907  }
908  }
909  } // end of mixed events loop
910  } // end of checking to see whether pool is filled correctly
911  pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
912  if(fDebug > 0 ) pool->PrintInfo();
913 }
914 //_____________________________________________________________________________
916 {
917  // terminate
918  if(fDebug > 0) cout << " *** Terminate() *** " << endl;
919 }
920 //______________________________________________________________________________
922 {
923  // Construct Flow Candidate Track from two selected candidates
924  if(fDebug > 1 ) cout << " *** MakeTrack() *** " << endl;
925  // if requested, check rapidity at this point
926  Bool_t overwrite = kTRUE;
927  AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
928  if (!sTrack) {
929  sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
930  overwrite = kFALSE;
931  }
932  else sTrack->ClearMe();
933  sTrack->SetMass(mass);
934  sTrack->SetPt(pt);
935  sTrack->SetPhi(phi);
936  sTrack->SetEta(eta);
937  sTrack->SetForPOISelection(kTRUE);
938  sTrack->SetForRPSelection(kFALSE);
939  if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
940  else fCandidates->AddLast(sTrack);
941  return;
942 }
943 //_____________________________________________________________________________
945 {
946  // Fill QA histos for MC analysis
947  TClonesArray *arrayMC = 0;
948  if(fDebug > 0) cout << " -> Switching to MC mode <- " << endl;
949  // fill array with mc tracks
950  arrayMC = (TClonesArray*) fAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
951  if (!arrayMC) AliFatal("Error: MC particles branch not found!\n");
952  for (Int_t iTracks = 0; iTracks < fAOD->GetNumberOfTracks(); iTracks++) {
953  AliAODTrack* track = dynamic_cast<AliAODTrack*>(fAOD->GetTrack(iTracks));
954  if(!track) AliFatal("Not a standard AOD");
955  if (fDebug>1) cout << " Received MC kaon " << endl;
956  Double_t b[2] = { -99., -99.};
957  Double_t bCov[3] = { -99., -99., -99.};
958  AliAODTrack copy(*track);
959  if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
960  // find corresponding mc particle
961  AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
962  if (!partMC) {
963  if(fDebug > 1) cout << "Cannot get MC particle" << endl;
964  continue;
965  }
966  // Check if it is primary, secondary from material or secondary from weak decay
967  Bool_t isPrimary = partMC->IsPhysicalPrimary();
968  Bool_t isSecondaryMaterial = kFALSE;
969  Bool_t isSecondaryWeak = kFALSE;
970  if (!isPrimary) {
971  Int_t mfl = -999, codemoth = -999;
972  Int_t indexMoth = partMC->GetMother();
973  if (indexMoth >= 0) { // is not fake
974  AliAODMCParticle* moth = (AliAODMCParticle*) arrayMC->At(indexMoth);
975  codemoth = TMath::Abs(moth->GetPdgCode());
976  mfl = Int_t(codemoth / TMath::Power(10, Int_t(TMath::Log10(codemoth))));
977  }
978  if (mfl == 3) isSecondaryWeak = kTRUE;
979  else isSecondaryMaterial = kTRUE;
980  }
981  if (isPrimary) {
982  fDCAPrim->Fill(track->Pt(), b[0]);
983  fDCAXYQA->Fill(b[0]);
984  fDCAZQA->Fill(b[1]);
985  }
986  if (isSecondaryWeak) fDCASecondaryWeak->Fill(track->Pt(), b[0]);
987  if (isSecondaryMaterial) fDCAMaterial->Fill(track->Pt(), b[0]);
988  }
989 }
990 
991 
992 
993 //_____________________________________________________________________________
994 Double_t AliAnalysisTaskCascadeTester::GetWDist(const AliVVertex* v0, const AliVVertex* v1)
995 {
996 
997  // calculate sqrt of weighted distance to other vertex
998  if (!v0 || !v1) {
999  printf("One of vertices is not valid\n");
1000  return 0;
1001  }
1002  static TMatrixDSym vVb(3);
1003  double dist = -1;
1004  double dx = v0->GetX()-v1->GetX();
1005  double dy = v0->GetY()-v1->GetY();
1006  double dz = v0->GetZ()-v1->GetZ();
1007  double cov0[6],cov1[6];
1008  v0->GetCovarianceMatrix(cov0);
1009  v1->GetCovarianceMatrix(cov1);
1010  //
1011  // fQxavsV0[0] fQxnmV0A
1012  // fQyavsV0[0] fQynmV0A
1013  // fQxavsV0[1] fQxnsV0A
1014  // fQyavsV0[1] fQynsV0A
1015  // fQxcvsV0[0] fQxnmV0C
1016  // fQycvsV0[0] fQynmV0C
1017  // fQxcvsV0[1] fQxnsV0C
1018  // fQycvsV0[1] fQynsV0C vVb(0,0) = cov0[0]+cov1[0];
1019  vVb(1,1) = cov0[2]+cov1[2];
1020  vVb(2,2) = cov0[5]+cov1[5];
1021  vVb(1,0) = vVb(0,1) = cov0[1]+cov1[1];
1022  vVb(0,2) = vVb(1,2) = vVb(2,0) = vVb(2,1) = 0.;
1023  vVb.InvertFast();
1024  if (!vVb.IsValid()) {printf("Singular Matrix\n"); return dist;}
1025  dist = vVb(0,0)*dx*dx + vVb(1,1)*dy*dy + vVb(2,2)*dz*dz
1026  + 2*vVb(0,1)*dx*dy + 2*vVb(0,2)*dx*dz + 2*vVb(1,2)*dy*dz;
1027  return dist>0 ? TMath::Sqrt(dist) : -1;
1028 
1029 }
1030 
1031 
1032 
1033 
1034 //_____________________________________________________________________________
1036 {
1037  // check for multi-vertexer pile-up
1038  //
1039  const int kMinPlpContrib = 5;
1040  const double kMaxPlpChi2 = 5.0;
1041  const double kMinWDist = 15;
1042  //
1043  const AliVVertex* vtPrm = 0;
1044  const AliVVertex* vtPlp = 0;
1045  int nPlp = 0;
1046  //
1047  if ( !(nPlp=aod->GetNumberOfPileupVerticesTracks()) ) return kFALSE;
1048  vtPrm = aod->GetPrimaryVertex();
1049  if (vtPrm == aod->GetPrimaryVertexSPD()) return kTRUE; // there are pile-up vertices but no primary
1050 
1051  //int bcPrim = vtPrm->GetBC();
1052  //
1053  for (int ipl=0;ipl<nPlp;ipl++) {
1054  vtPlp = (const AliVVertex*)aod->GetPileupVertexTracks(ipl);
1055  //
1056  if (vtPlp->GetNContributors() < kMinPlpContrib) continue;
1057  if (vtPlp->GetChi2perNDF() > kMaxPlpChi2) continue;
1058  // int bcPlp = vtPlp->GetBC();
1059  // if (bcPlp!=AliVTrack::kTOFBCNA && TMath::Abs(bcPlp-bcPrim)>2) return kTRUE; // pile-up from other BC
1060  //
1061  double wDst = GetWDist(vtPrm,vtPlp);
1062  if (wDst<kMinWDist) continue;
1063  //
1064  return kTRUE; // pile-up: well separated vertices
1065  }
1066  //
1067  return kFALSE;
1068  //
1069 }
1070 
1071 //_____________________________________________________________________________
1072 //---------------------------------------------------
1074 {
1075 
1076  Short_t kPID = 0;
1077 
1078  if((nSk == nSpi) && (nSk == nSp))
1079  return kPID;
1080  if((nSk < nSpi) && (nSk < nSp) && (nSk < fPIDConfig[0]))
1081  kPID = 2;
1082 
1083  if((nSpi < nSk) && (nSpi < nSp) && (nSpi < fPIDConfig[0]))
1084  kPID = 1;
1085 
1086  if((nSp < nSk) && (nSp < nSpi) && (nSp < fPIDConfig[0]))
1087  kPID = 3;
1088  return kPID;
1089 
1090 }
1091 
1092 
1093 //_____________________________________________________________________________
1095 {
1096 
1097  Bool_t kDC = kFALSE;
1098 
1099  if (nSk < fPIDConfig[0] && minNSigma != 2)
1100  kDC = kTRUE;
1101 
1102  return kDC;
1103 
1104 }
1105 
const Color_t cc[]
Definition: DrawKs.C:1
double Double_t
Definition: External.C:58
void SetEta(Double_t eta)
TH1F * fInvMPP[18]
like-sign kaon pairs
Definition: External.C:236
TH2F * fNOPIDTOF
QA histogram of TPC response of kaons.
Bool_t fSkipEventSelection
profiles for vzero vn(minv)
TH1F * InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
TH2F * fPIDk
QA histogram of TPC response of all charged particles.
void PlotMultiplcities(const T *event) const
Int_t GetNumberOfRPs() const
TH2F * fMultCorAfterCuts
QA histogram of p_t distribution of negative kaons.
TH1F * fVZEROA
QA plot of eta distribution of tracks used for event plane estimation.
Double_t mass
TProfile * fV0Data[18][2]
subevent resolution info for v2
TH1F * fTPCM
QA plot vzeroc multiplicity (all tracks in event)
TH1F * fVZEROC
QA plot vzeroa multiplicity (all tracks in event)
TH1F * fPtKP
QA histogram of p_t distribution of negative particles.
TCanvas * c
Definition: TestFitELoss.C:172
Double_t PhiPt(const T *track_1, const T *track_2) const
TH1F * fInvMNP[18]
QA histo of TOF response kaons.
void SetMass(Double_t mass)
TH2F * BookPIDHistogram(const char *name, Bool_t TPC)
virtual Bool_t IsSelected(TObject *obj, Int_t id=-666)
void SetParamType(trackParameterType paramType)
TH2F * fDCAAll
QA plot TPC multiplicity (tracks used for event plane estimation)
void ReconstructionWithEventMixing(TObjArray *MixingCandidates) const
void SetReferenceMultiplicity(Int_t m)
Double_t Mass() const
TH1F * fPtN
QA histogram of p_t distribution of positive particles.
AliESDv0 * fV0
placeholder for the current kink
Bool_t InRPSelection() const
AliFlowBayesianPID * fBayesianResponse
AliFlowTrack * GetTrack(Int_t i)
AliEventPoolManager * fPoolManager
AOD oject.
UShort_t T(UShort_t m, UShort_t t)
Definition: RingBits.C:60
void SetForRPSelection(Bool_t b=kTRUE)
void SetPtRange(Float_t r1, Float_t r2)
Int_t GetNDaughters(void) const
void SetNumberOfRPs(Int_t nr)
void Fill(AliFlowTrackCuts *rpCuts, AliFlowTrackCuts *poiCuts)
TH1F * fPtSpectra[18]
like-sign kaon pairs
int Int_t
Definition: External.C:63
Double_t InvariantMass(const T *track1, const T *track2) const
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
TH2F * fDCASecondaryWeak
dca of primaries (mc) or kaons (data)
Bool_t CheckCandidateEtaPtCut(const T *track1, const T *track2) const
void MakeTrack(Double_t, Double_t, Double_t, Double_t) const
Bool_t plpMV(const AliAODEvent *aod)
static AliFlowCommonConstants * GetMaster()
Double_t fCentralityMin
QA profile of centralty vs multiplicity.
TH1F * fPt
QA plot of azimuthal distribution of tracks used for event plane estimation.
Bool_t GetDoubleCountingK(Double_t nSk, Short_t minNSigma) const
TH1F * fEta
QA plot of p_t sectrum of tracks used for event plane estimation.
TProfile * fSubEventDPhiv2
dca material (mc) all (data)
AliPIDResponse * fPIDResponse
void SetPhi(Double_t phi)
short Short_t
Definition: External.C:23
void PtSelector(Int_t _track_type, const T *track_1, const T *track_2) const
void DefineDeadZone(Double_t etaMin, Double_t etaMax, Double_t phiMin, Double_t phiMax)
void SetPt(Double_t pt)
Int_t GetIDDaughter(Int_t value) const
Short_t FindMinNSigma(Double_t nSpi, Double_t nSk, Double_t nSp) const
void SetEvent(AliVEvent *event, AliMCEvent *mcEvent=NULL)
Float_t nEvents[nProd]
TH1F * fPtKN
QA histogram of p_t distribution of positive kaons.
TH2F * fNOPID
QA histogram of events that do not pass centrality cut.
TH2F * fMultvsCentr
QA profile global and tpc multiplicity after outlier cut.
Double_t GetWDist(const AliVVertex *v0, const AliVVertex *v1)
void SetForPOISelection(Bool_t b=kTRUE)
TList * fOutputList
event pool manager
virtual Int_t Charge() const
const char Option_t
Definition: External.C:48
virtual void UserExec(Option_t *option)
bool Bool_t
Definition: External.C:53
TObjArray * fCandidates
PID response object.
TH2F * fPIDTOF
QA histo of TOF repsonse charged particles.
TH1F * fInvMNN[18]
unlike sign kaon pairs
void InsertTrack(AliFlowTrack *)
virtual void ClearFast()
ClassImp(AliAnalysisTaskCascadeTester) AliAnalysisTaskCascadeTester
AliFlowEvent * fFlowEvent
pid response object
TH1F * fDCAXYQA
qa dca of all charged particles
void SetEtaRange(Float_t r1, Float_t r2)
Int_t NumberOfTracks() const