]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGCF/FLOW/Attic/AliAnalysisTwoParticleResonanceFlowTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / FLOW / Attic / AliAnalysisTwoParticleResonanceFlowTask.cxx
CommitLineData
fbd0da9c 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// AliAnalysisTwoParticleResonanceFlowTask:
16// author: Redmer Alexander Bertens (rbertens@cern.ch)
17// You Zhou (you.zhou@cern.ch)
18// analyis task for
19// AliResonanceFlowHelperTrack provides a lightweight helper track for reconstruction
20// This analysis DOES NOT support ESDs
21// Version: v2012.9.7
22
23
fbd0da9c 24#include "TH1F.h"
25#include "TH2F.h"
fbd0da9c 26#include "TMath.h"
9e437f94 27#include "TObject.h"
fbd0da9c 28#include "TObjArray.h"
fbd0da9c 29#include "AliAnalysisTaskSE.h"
30#include "AliAnalysisManager.h"
31#include "AliAODEvent.h"
32#include "AliAODInputHandler.h"
33#include "AliCentrality.h"
fbd0da9c 34#include "AliAnalysisTwoParticleResonanceFlowTask.h"
35#include "AliFlowBayesianPID.h"
36#include "AliStack.h"
37#include "AliMCEvent.h"
38#include "TProfile.h"
9e437f94 39#include "TDirectoryFile.h"
fbd0da9c 40#include "AliFlowCandidateTrack.h"
41#include "AliFlowTrackCuts.h"
f5cd76c4 42#include "AliFlowEventCuts.h"
fbd0da9c 43#include "AliFlowEventSimple.h"
44#include "AliFlowTrackSimple.h"
45#include "AliFlowCommonConstants.h"
46#include "AliFlowEvent.h"
47#include "TVector3.h"
48#include "AliAODVZERO.h"
49#include "AliPID.h"
50#include "AliPIDResponse.h"
51#include "AliAODMCParticle.h"
52#include "AliAnalysisTaskVnV0.h"
53#include "AliEventPoolManager.h"
54
55class AliFlowTrackCuts;
56
57using std::cout;
58using std::endl;
59
60ClassImp(AliAnalysisTwoParticleResonanceFlowTask)
61
62AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask() : AliAnalysisTaskSE(),
f5cd76c4 63 fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fCutsEvent(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fMultCorAfterCuts(0), fMultvsCentr(0), fPhi(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), fAnalysisSummary(0)
fbd0da9c 64{
65 // Default constructor
66 for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
67 for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
68 for(Int_t i(0); i < 20; i++) {
69 fVertexMixingBins[i] = 0;
70 fCentralityMixingBins[i] = 0;
71 }
72 fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
73 for(Int_t i(0); i < 18; i++) {
74 for(Int_t j(0); j < 2; j++) {
75 fV0Data[i][j] = 0;
76 for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
77 for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
78 }
79 fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
80 }
70adca5c 81 for(Int_t i(0); i < 12; i++) fAddTaskMacroSummary[i] = 0.;
fbd0da9c 82}
83//_____________________________________________________________________________
f5cd76c4 84AliAnalysisTwoParticleResonanceFlowTask::AliAnalysisTwoParticleResonanceFlowTask(const char *name) : AliAnalysisTaskSE(name), fSpeciesA(0), fSpeciesB(0), fChargeA(0), fChargeB(0), fMassA(0), fMassB(0), fMinPtA(0), fMaxPtA(0), fMinPtB(0), fMaxPtB(0), fIsMC(0), fEventMixing(0), fPhiMinusPsiMethod(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fCutsEvent(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fNdPhiBins(18), fCentrality(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fPIDp(0), fPtP(0), fPtN(0), fPtSpeciesA(0), fPtSpeciesB(0), fMultCorAfterCuts(0), fMultvsCentr(0), fPhi(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), fAnalysisSummary(0)
fbd0da9c 85{
86 // Constructor
87 for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
88 for(Int_t i(0); i < 5; i++) fDCAConfig[i] = 0.;
89 for(Int_t i(0); i < 20; i++) {
90 fVertexMixingBins[i] = 0;
91 fCentralityMixingBins[i] = 0;
92 }
93 fMixingParameters[0] = 1000; fMixingParameters[1] = 50000; fMixingParameters[2] = 5;
94 for(Int_t i(0); i < 18; i++) {
95 for(Int_t j(0); j < 2; j++) {
96 fV0Data[i][j] = 0;
97 for(Int_t k(0); k < 15; k++) fPhiMinusPsiDataContainer[k][i][j] = 0x0;
98 for(Int_t k(0); k < 15; k++) fPhiMinusPsiBackgroundContainer[k][i][j] = 0x0;
99 }
100 fResonanceSignal[i] = 0; fResonanceBackground[i] = 0; fPtSpectra[i] = 0; fPtBins[i] = 0.; fdPhiBins[i] = 0.;
101 }
70adca5c 102 for(Int_t i(0); i < 12; i++) fAddTaskMacroSummary[i] = 0.;
fbd0da9c 103 DefineInput(0, TChain::Class());
104 DefineOutput(1, TList::Class());
105 DefineOutput(2, AliFlowEventSimple::Class());
fbd0da9c 106}
107//_____________________________________________________________________________
108AliAnalysisTwoParticleResonanceFlowTask::~AliAnalysisTwoParticleResonanceFlowTask()
109{
110 // Destructor
111 if (fNullCuts) delete fNullCuts;
112 if (fOutputList) delete fOutputList;
113 if (fCandidates) delete fCandidates;
114 if (fFlowEvent) delete fFlowEvent;
115 if (fBayesianResponse) delete fBayesianResponse;
116 if (fEventMixing) delete fPoolManager;
fbd0da9c 117}
118//_____________________________________________________________________________
119void AliAnalysisTwoParticleResonanceFlowTask::ForceExit(Int_t type, const char* message)
120{
121 // force exit
122 if(type==0) cout << " --> Illegal options in AddTask <-- " << endl;
123 if(type==1) cout << " --> Unknown error <-- " << endl;
124 AliFatal(message);
125}
126//_____________________________________________________________________________
127TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookHistogram(const char* name)
128{
129 // Return a pointer to a TH1 with predefined binning
fbd0da9c 130 TH1F *hist = new TH1F(name, Form("M_{INV} (%s)", name), 2*fMassBins, fMinMass, fMaxMass);
131 hist->GetXaxis()->SetTitle("M_{INV} (GeV / c^{2})");
132 hist->GetYaxis()->SetTitle("No. of pairs");
133 hist->SetMarkerStyle(kFullCircle);
134 hist->Sumw2();
135 fOutputList->Add(hist);
136 return hist;
137}
138//_____________________________________________________________________________
139TH2F* AliAnalysisTwoParticleResonanceFlowTask::BookPIDHistogram(const char* name, Bool_t TPC)
140{
141 // Return a pointer to a TH2 with predefined binning
fbd0da9c 142 TH2F *hist = 0x0;
143 if(TPC) {
144 hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1000);
145 hist->GetYaxis()->SetTitle("dE/dx (a.u.)");
146 }
147 if(!TPC) {
148 hist = new TH2F(name, Form("PID (%s)", name), 100, 0, 5, 100, 0, 1.5);
149 hist->GetYaxis()->SetTitle("#beta");
150 }
151 hist->GetXaxis()->SetTitle("P (GeV / c)");
152 fOutputList->Add(hist);
153 return hist;
154}
155//_____________________________________________________________________________
156TH1F* AliAnalysisTwoParticleResonanceFlowTask::InitPtSpectraHistograms(Float_t nmin, Float_t nmax)
157{
158 // intialize p_t histograms for each p_t bin
fbd0da9c 159 TH1F* hist = new TH1F(Form("%4.2f p_{t} %4.2f", nmin, nmax), Form("%f p_{t} %f", nmin, nmax), 2*fMassBins, nmin, nmax);
160 hist->GetXaxis()->SetTitle("p_{T} GeV / c");
161 fOutputList->Add(hist);
162 return hist;
163}
164//_____________________________________________________________________________
165TH1F* AliAnalysisTwoParticleResonanceFlowTask::BookPtHistogram(const char* name)
166{
167 // Return a pointer to a p_T spectrum histogram
fbd0da9c 168 TH1F* ratio = new TH1F(name, name, 100, 0, 7);
169 ratio->GetXaxis()->SetTitle("p_{T} ( GeV / c^{2} )");
170 ratio->GetYaxis()->SetTitle("No. of events");
171 ratio->Sumw2();
172 fOutputList->Add(ratio);
173 return ratio;
174}
175//_____________________________________________________________________________
176Bool_t AliAnalysisTwoParticleResonanceFlowTask::InitializeAnalysis()
177{
178 // set up the analysis
179 fNullCuts = new AliFlowTrackCuts("null_cuts");
180 fBayesianResponse = new AliFlowBayesianPID();
181 Float_t t(0);
182 for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
183 if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
184 if(fNPtBins > 18) AliFatal("Invalid number of pt bins initialied ( > 18 ) -- terminating analysis !!!");
185 AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
186 cc->SetNbinsQ(500); cc->SetNbinsPhi(180); cc->SetNbinsMult(10000);
187 cc->SetQMin(0.0); cc->SetPhiMin(0.0); cc->SetMultMin(0);
188 cc->SetQMax(3.0); cc->SetPhiMax(TMath::TwoPi()); cc->SetMultMax(10000);
189 cc->SetNbinsMass(fMassBins); cc->SetNbinsEta(200); (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
190 cc->SetMassMin(fMinMass); cc->SetEtaMin(-5.0); cc->SetPtMin(0);
191 cc->SetMassMax(fMaxMass); cc->SetEtaMax(+5.0); (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
192 AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
193 if (man) {
194 AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
195 if (inputHandler) fPIDResponse = inputHandler->GetPIDResponse();
196 }
197 fMassA = fMassA*fMassA; //we always use mass squared in calculations, so doing it once in the init reduces cpu load
198 fMassB = fMassB*fMassB;
199 return kTRUE;
200}
201//_____________________________________________________________________________
202void AliAnalysisTwoParticleResonanceFlowTask::UserCreateOutputObjects()
203{
204 // Create user defined output objects
205 if(!InitializeAnalysis()) ForceExit(0, " > Couldn't initialize the analysis !!! <");
fbd0da9c 206 // Create all output objects and store them to a list
207 fOutputList = new TList();
208 fOutputList->SetOwner(kTRUE);
209 if(fQA) {
210 fEventStats = new TH1F("fHistStats", "Event Statistics", 18, -.25, 4.25);
211 fEventStats->GetXaxis()->SetTitle("No. of events");
212 fEventStats->GetYaxis()->SetTitle("Statistics");
213 fOutputList->Add(fEventStats);
214 }
215 fCentralityPass = new TH1F("fCentralityPass", "Centrality Pass", 101, -1, 100);
216 fOutputList->Add(fCentralityPass);
217 if(fQA) {
218 fCentralityNoPass = new TH1F("fCentralityNoPass", "Centrality No Pass", 101, -1, 100);
219 fOutputList->Add(fCentralityNoPass);
220 fNOPID = BookPIDHistogram("TPC signal, all particles", kTRUE);
8da4b0c9 221 fPIDk = BookPIDHistogram("TPC signal, Species A", kTRUE);
222 fPIDp = BookPIDHistogram("TPC signal, Species B", kTRUE);
fbd0da9c 223 }
224 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
225 if(!fPhiMinusPsiMethod) {
226 fResonanceSignal[ptbin] = BookHistogram(Form("NP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
227 fResonanceBackground[ptbin] = BookHistogram(Form("PP, %4.2f < p_{T} < %4.2f GeV", fPtBins[ptbin], fPtBins[ptbin+1]));
228 }
70adca5c 229 else if(fPhiMinusPsiMethod) {
fbd0da9c 230 for(Int_t i(0); i<fNdPhiBins; i++) {// for all delta phi bins
231 fPhiMinusPsiDataContainer[i][ptbin][0] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
232 fPhiMinusPsiDataContainer[i][ptbin][1] = BookHistogram(Form("%4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
233 fPhiMinusPsiBackgroundContainer[i][ptbin][0] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{A}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
234 fPhiMinusPsiBackgroundContainer[i][ptbin][1] = BookHistogram(Form("BG, %4.2f < p_{T} < %4.2f GeV, %4.2f < (#phi - #Psi_{C}) < %4.2f", fPtBins[ptbin], fPtBins[ptbin+1], fdPhiBins[i], fdPhiBins[i+1]));
235 }
236 }
237 }
238 if(fQA) {
239 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) fPtSpectra[ptbin] = InitPtSpectraHistograms(fPtBins[ptbin], fPtBins[ptbin+1]);
240 fPtP = BookPtHistogram("i^{+}");
241 fPtN = BookPtHistogram("i^{-}");
8da4b0c9 242 fPtSpeciesA = BookPtHistogram("p_{T} spectrum Species A");
243 fPtSpeciesB = BookPtHistogram("p_{T} spectrum Species B");
fbd0da9c 244 fPhi = new TH1F("fPhi", "#phi distribution", 100, -.5, 7);
245 fOutputList->Add(fPhi);
fbd0da9c 246 fEta = new TH1F("fEta", "#eta distribution", 100, -1.1, 1.1);
247 fOutputList->Add(fEta);
248 fVZEROA = new TH1F("fVZEROA", "VZERO A Multiplicity", 1000, 0, 10000);
249 fOutputList->Add(fVZEROA);
250 fVZEROC = new TH1F("fVZEROC", "VZERO C Multiplicity", 1000, 0, 10000);
251 fOutputList->Add(fVZEROC);
252 fTPCM = new TH1F("fTPCM", "TPC multiplicity", 1000, 0, 10000);
253 fOutputList->Add(fTPCM);
254 fDCAXYQA = new TH1F("fDCAXYQA", "fDCAXYQA", 1000, -5, 5);
255 fOutputList->Add(fDCAXYQA);
256 fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
257 fOutputList->Add(fDCAZQA);
089d87ea 258 if(fCentralityCut2010 || fCentralityCut2011) {
90555064 259 fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
260 fOutputList->Add(fMultCorAfterCuts);
261 fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
262 fOutputList->Add(fMultvsCentr);
263 }
fbd0da9c 264 }
265 if(fIsMC || fQA) {
266 fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
267 fOutputList->Add(fDCAAll);
268 fDCAPrim = new TH2F("fDCAprim","fDCAprim", 1000, 0, 10, 1000, -5, 5);
269 fOutputList->Add(fDCAPrim);
270 fDCASecondaryWeak = new TH2F("fDCASecondaryWeak","fDCASecondaryWeak", 1000, 0, 10, 1000, -5, 5);
271 fOutputList->Add(fDCASecondaryWeak);
272 fDCAMaterial = new TH2F("fDCAMaterial","fDCAMaterial", 1000, 0, 10, 1000, -5, 5);
273 fOutputList->Add(fDCAMaterial);
274 }
275 if(fV0||fPhiMinusPsiMethod) {
276 fSubEventDPhiv2 = new TProfile("fSubEventDPhiv2", "fSubEventDPhiv2", 3, 0, 3);
277 fSubEventDPhiv2->GetXaxis()->SetBinLabel(1, "<cos(2(#Psi_{a} - #Psi_{b}))>");
278 fSubEventDPhiv2->GetXaxis()->SetBinLabel(2, "<cos(2(#Psi_{a} - #Psi_{c}>))");
279 fSubEventDPhiv2->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{b} - #Psi_{c}>))");
280 fOutputList->Add(fSubEventDPhiv2);
281 const char* V0[] = {"V0A", "V0C"};
282 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++)
283 for(Int_t iV0(0); iV0 < 2; iV0++) {
284 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);
285 fOutputList->Add(fV0Data[ptbin][iV0]);
286 }
287 }
8da4b0c9 288 TString name = "fAnalysisSummary_";
289 name+=this->GetName();
290 fAnalysisSummary = new TH1F(name.Data(), name.Data(), 34, 0, 34);
fbd0da9c 291 fAnalysisSummary->GetXaxis()->SetBinLabel(1, "fIsMC");
292 fAnalysisSummary->GetXaxis()->SetBinLabel(2, "fEventMixing");
293 fAnalysisSummary->GetXaxis()->SetBinLabel(3, "fAODAnalysis");
8da4b0c9 294 if(fPIDConfig[6] > 0) fAnalysisSummary->GetXaxis()->SetBinLabel(4, "bayesian purity");
295 else fAnalysisSummary->GetXaxis()->SetBinLabel(4, "TPC TOF PID");
fbd0da9c 296 fAnalysisSummary->GetXaxis()->SetBinLabel(5, "dca mode");
297 fAnalysisSummary->GetXaxis()->SetBinLabel(6, "fVertex");
298 fAnalysisSummary->GetXaxis()->SetBinLabel(7, "fCentralityMin");
299 fAnalysisSummary->GetXaxis()->SetBinLabel(8, "fCentralityMax");
8da4b0c9 300 fAnalysisSummary->GetXaxis()->SetBinLabel(9, "centrality");
1164b0ff 301 fAnalysisSummary->GetXaxis()->SetBinLabel(10, fkCentralityMethodA);
fbd0da9c 302 fAnalysisSummary->GetXaxis()->SetBinLabel(11, "fApplyDeltaDipCut");
303 fAnalysisSummary->GetXaxis()->SetBinLabel(12, "POI filterbit");
304 fAnalysisSummary->GetXaxis()->SetBinLabel(13, "RP filterbit");
305 fAnalysisSummary->GetXaxis()->SetBinLabel(14, "PhiMinusPsiMethod");
306 fAnalysisSummary->GetXaxis()->SetBinLabel(15, "SP");
307 fAnalysisSummary->GetXaxis()->SetBinLabel(16, "SPSUB");
308 fAnalysisSummary->GetXaxis()->SetBinLabel(17, "QC");
309 fAnalysisSummary->GetXaxis()->SetBinLabel(18, "EP");
310 fAnalysisSummary->GetXaxis()->SetBinLabel(19, "EP3sub");
311 fAnalysisSummary->GetXaxis()->SetBinLabel(20, "V0_SP");
312 fAnalysisSummary->GetXaxis()->SetBinLabel(21, "eta gap");
313 fAnalysisSummary->GetXaxis()->SetBinLabel(22, "harm");
314 fAnalysisSummary->GetXaxis()->SetBinLabel(23, "highPtMode");
315 fAnalysisSummary->GetXaxis()->SetBinLabel(24, "deltaDipAngle");
8da4b0c9 316 fAnalysisSummary->GetXaxis()->SetBinLabel(26, "SpeciesA");
317 fAnalysisSummary->GetXaxis()->SetBinLabel(27, "chargeA");
318 fAnalysisSummary->GetXaxis()->SetBinLabel(28, "fMinPtA");
319 fAnalysisSummary->GetXaxis()->SetBinLabel(29, "fMaxPtA");
320 fAnalysisSummary->GetXaxis()->SetBinLabel(30, "SpeciesB");
321 fAnalysisSummary->GetXaxis()->SetBinLabel(31, "chargeB");
322 fAnalysisSummary->GetXaxis()->SetBinLabel(32, "fMinPtB");
323 fAnalysisSummary->GetXaxis()->SetBinLabel(33, "fMaxPtB");
324 fAnalysisSummary->GetXaxis()->SetBinLabel(34, "iterator");
fbd0da9c 325 fOutputList->Add(fAnalysisSummary);
326 fAnalysisSummary->SetBinContent(1, (Float_t)fIsMC);
327 fAnalysisSummary->SetBinContent(2, (Float_t)fEventMixing);
328 fAnalysisSummary->SetBinContent(3, (Float_t)1);
8da4b0c9 329 if(fPIDConfig[6] > 0) fAnalysisSummary->SetBinContent(4, (Float_t)fPIDConfig[6]);
fbd0da9c 330 fAnalysisSummary->SetBinContent(5, (Float_t)fDCAConfig[0]);
331 fAnalysisSummary->SetBinContent(6, (Float_t)fVertexRange);
332 fAnalysisSummary->SetBinContent(7, (Float_t)fCentralityMin);
333 fAnalysisSummary->SetBinContent(8, (Float_t)fCentralityMax);
fbd0da9c 334 fAnalysisSummary->SetBinContent(11, (Float_t)fApplyDeltaDipCut);
335 fAnalysisSummary->SetBinContent(12, (Float_t)fPOICuts->GetAODFilterBit());
336 fAnalysisSummary->SetBinContent(13, (Float_t)fCutsRP->GetAODFilterBit());
70adca5c 337 for(Int_t i(0); i<12;i++) fAnalysisSummary->SetBinContent(14+i, fAddTaskMacroSummary[i]);
8da4b0c9 338 fAnalysisSummary->SetBinContent(26, (Float_t)fSpeciesA);
339 fAnalysisSummary->SetBinContent(27, (Float_t)fChargeA);
340 fAnalysisSummary->SetBinContent(28, (Float_t)fMinPtA);
341 fAnalysisSummary->SetBinContent(29, (Float_t)fMaxPtA);
342 fAnalysisSummary->SetBinContent(30, (Float_t)fSpeciesB);
343 fAnalysisSummary->SetBinContent(31, (Float_t)fChargeB);
344 fAnalysisSummary->SetBinContent(32, (Float_t)fMinPtB);
345 fAnalysisSummary->SetBinContent(33, (Float_t)fMaxPtB);
346 fAnalysisSummary->SetBinContent(34, (Float_t)1);
fbd0da9c 347 PostData(1, fOutputList);
348 // create candidate array
349 fCandidates = new TObjArray(1000);
350 fCandidates->SetOwner(kTRUE);
351 // create and post flowevent
352 fFlowEvent = new AliFlowEvent(10000);
353 PostData(2, fFlowEvent);
354 if(fEventMixing) fPoolManager = InitializeEventMixing();
355}
356//_____________________________________________________________________________
357AliEventPoolManager* AliAnalysisTwoParticleResonanceFlowTask::InitializeEventMixing()
358{
359 // initialize event mixing
fbd0da9c 360 Int_t _c(0), _v(0);
361 for(Int_t i(0); i < 19; i++) {
362 if (fCentralityMixingBins[i+1] < fCentralityMixingBins[i]) { _c = i; break; }
363 else _c = 19;
364 }
365 for(Int_t i(0); i < 19; i++) {
366 if (fVertexMixingBins[i+1] < fVertexMixingBins[i]) { _v = i; break; }
367 else _v = 19;
368 }
63c84ccc 369 Double_t centralityBins[_c];
370 Double_t vertexBins[_v];
90555064 371 for(Int_t i(0); i < _c + 1; i++) centralityBins[i] = fCentralityMixingBins[i];
372 for(Int_t i(0); i < _v + 1; i++) vertexBins[i] = fVertexMixingBins[i];
373 return new AliEventPoolManager(fMixingParameters[0], fMixingParameters[1], _c, (Double_t*)centralityBins, _v, (Double_t*)vertexBins);
fbd0da9c 374}
375//_____________________________________________________________________________
376template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::InvariantMass(const T* track1, const T* track2) const
377{
378 // Return the invariant mass of two tracks
fbd0da9c 379 if ((!track2) || (!track1)) return 0.;
380 Float_t pxs = TMath::Power((track1->Px() + track2->Px()), 2);
381 Float_t pys = TMath::Power((track1->Py() + track2->Py()), 2);
382 Float_t pzs = TMath::Power((track1->Pz() + track2->Pz()), 2);
383 Float_t e1 = TMath::Sqrt(track1->P() * track1->P() + track1->Mass());
384 Float_t e2 = TMath::Sqrt(track2->P() * track2->P() + track2->Mass());
385 Float_t es = TMath::Power((e1 + e2), 2);
386 if ((es - (pxs + pys + pzs)) < 0) return 0.;
387 return TMath::Sqrt((es - (pxs + pys + pzs)));
388}
389//_____________________________________________________________________________
390template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::DeltaDipAngle(const T* track1, const T* track2) const
391{
392 // Calculate the delta dip angle between two particles
fbd0da9c 393 if (track1->P()*track2->P() == 0) return 999;
394 return TMath::ACos(((track1->Pt() * track2->Pt()) + (track1->Pz() * track2->Pz())) / (track1->P() * track2->P()));
395}
396//_____________________________________________________________________________
397template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckDeltaDipAngle(const T* track1, const T* track2) const
398{
399 // Check if pair passes delta dip angle cut within 0 < p_t < fDeltaDipPt
fbd0da9c 400 if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PairPt(track1, track2) < fDeltaDipPt)) return kFALSE;
401 return kTRUE;
402}
403//_____________________________________________________________________________
404template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
405{
406 // Check if pair passes eta and pt cut
fbd0da9c 407 if (fCandidateMinPt > PairPt(track1, track2) || fCandidateMaxPt < PairPt(track1, track2)) return kFALSE;
408 TVector3 a(track1->Px(), track1->Py(), track1->Pz());
409 TVector3 b(track2->Px(), track2->Py(), track2->Pz());
410 TVector3 c = a + b;
411 if (fCandidateMinEta > c.Eta() || fCandidateMaxEta < c.Eta()) return kFALSE;
412 return kTRUE;
413}
414//_____________________________________________________________________________
fbd0da9c 415template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::PlotMultiplcities(const T* event) const
416{
417 // QA multiplicity plots
fbd0da9c 418 fVZEROA->Fill(event->GetVZEROData()->GetMTotV0A());
419 fVZEROC->Fill(event->GetVZEROData()->GetMTotV0C());
420 fTPCM->Fill(event->GetNumberOfTracks());
421}
422//_____________________________________________________________________________
fbd0da9c 423void AliAnalysisTwoParticleResonanceFlowTask::InitializeBayesianPID(AliAODEvent* event)
424{
425 // Initialize the Bayesian PID object for AOD
fbd0da9c 426 fBayesianResponse->SetDetResponse(event, fCentrality);
427}
428//_____________________________________________________________________________
429template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesTPCbayesianCut(T* track, Int_t species) const
430{
431 // Check if the particle passes the TPC TOF bayesian cut.
fbd0da9c 432 fBayesianResponse->ComputeProb(track);
433 if (!fBayesianResponse->GetCurrentMask(0)) return kFALSE; // return false if TPC has no response
434 Int_t contamination(0); //type of particle we don't want in the sample
435 if(species==3) contamination = 2; //we want to reject pions when selection kaons
436 if(species==2) contamination = 3; // we want to reject kaons when selecting pions
437 Float_t *probabilities = fBayesianResponse->GetProb();
8da4b0c9 438 return (probabilities[species] > fPIDConfig[6] && (probabilities[contamination] < probabilities[species]));// 3 is kaon, 2 is pion
fbd0da9c 439}
440//_____________________________________________________________________________
441Bool_t AliAnalysisTwoParticleResonanceFlowTask::PassesDCACut(AliAODTrack* track) const
442{
443 // check if track passes dca cut according to dca routine
444 // setup the routine as follows:
445 // fDCAConfig[0] < -1 no pt dependence
446 // fDCAConfig[0] = 0 do nothing
447 // fDCAConfig[0] > 1 pt dependent dca cut
646a809c 448 AliAODTrack copy(*track);
fbd0da9c 449 if(fIsMC) return kTRUE;
450 if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
451 if(fQA) {
452 Double_t b[2] = { -99., -99.};
453 Double_t bCov[3] = { -99., -99., -99.};
646a809c 454 if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
9ea73210 455 fDCAXYQA->Fill(b[0]);
456 fDCAZQA->Fill(b[1]);
457 fDCAPrim->Fill(track->Pt(), b[0]);
458 }
fbd0da9c 459 }
460 return kTRUE;
461 }
462 Double_t b[2] = { -99., -99.};
463 Double_t bCov[3] = { -99., -99., -99.};
646a809c 464 if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
fbd0da9c 465 if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
466 if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
467 if(fDCAConfig[0] > .9) {
468 if(fDCAConfig[4] < TMath::Abs(b[1])) return kFALSE;
469 Float_t denom = TMath::Power(track->Pt(), fDCAConfig[3]);
470 if( denom < 0.0000001 ) return kFALSE; // avoid division by zero
471 if( (fDCAConfig[1] + fDCAConfig[2] / denom) < TMath::Abs(b[0]) ) return kFALSE;
472 }
473 if(fQA) {
474 fDCAXYQA->Fill(b[0]);
475 fDCAZQA->Fill(b[1]);
476 fDCAPrim->Fill(track->Pt(), b[0]);
477 }
478 return kTRUE;
479}
480//_____________________________________________________________________________
8da4b0c9 481Bool_t AliAnalysisTwoParticleResonanceFlowTask::DoOwnPID(AliAODTrack* track, Int_t species) const
482{
483 // do custom pid based on tpc and tof response
484 Float_t nSigmasTPC = (species==fSpeciesA) ? fPIDConfig[0] : fPIDConfig[3]; // number that is used when only tpc is available (5)
485 Float_t nSigmasTPCwithTOF = (species==fSpeciesA) ? fPIDConfig[1] : fPIDConfig[4]; // number that is used for tpc when tof is available (5)
486 Float_t nSigmasTOF = (species==fSpeciesA) ? fPIDConfig[2] : fPIDConfig[5]; // number that is used for tof (if available) (3)
487 if ((track->GetStatus() & AliESDtrack::kTOFout)&&(track->GetStatus() & AliESDtrack::kTIME)) { // check if the track is matched by tof signal
488 if (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)species)) > nSigmasTPCwithTOF) return kFALSE; // reject with tpc
489 if (track->P() < 1.5) return (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)species)) <= nSigmasTOF);
490 return (TMath::Abs(fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)species)) <= (nSigmasTOF-.333*nSigmasTOF));
491 }
492 else { // switch to tpc pid if no tof hit is found
493 if(track->Pt() < .35) return (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)species)) <= nSigmasTPC);
494 else if(track->Pt() < .5) return (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)species)) < (nSigmasTPC-.4*nSigmasTPC));
495 return (TMath::Abs(fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)species)) <= (nSigmasTPC-.6*nSigmasTPC));
496 }
497 return kFALSE;
498}
499//_____________________________________________________________________________
fbd0da9c 500Bool_t AliAnalysisTwoParticleResonanceFlowTask::AcceptTrack(AliAODTrack* track, Int_t species) const
501{
8da4b0c9 502 if(fQA) TrackQA(track, species, kTRUE); // do qa before cuts
503 if(species==fSpeciesA&&((track->Pt()<fMinPtA)||(track->Pt()>fMaxPtA))) return kFALSE; // least expensive check first
504 if(species==fSpeciesB&&((track->Pt()<fMinPtB)||(track->Pt()>fMaxPtB))) return kFALSE;
505 if(!PassesDCACut(track)) return kFALSE;
506 if(fPIDConfig[6] < 0) {
507 if(DoOwnPID(track, species)) {
508 if(fQA) TrackQA(track, species, kFALSE);
509 return kTRUE;
510 }
511 return kFALSE;
512 }
513 else {
514 if (PassesTPCbayesianCut(track, species)) {
515 if(fQA) TrackQA(track, species, kFALSE);
516 return kTRUE;
517 }
fbd0da9c 518 }
519 return kFALSE;
520}
521//_____________________________________________________________________________
522template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PairPt(const T* track1, const T* track2, Bool_t phi) const
523{
524 // return p_t of track pair, select phi to return the azimuthal angle instead
525 TVector3 a(track1->Px(), track1->Py(), track1->Pz());
526 TVector3 b(track2->Px(), track2->Py(), track2->Pz());
527 TVector3 c = a + b;
528 if(phi) return c.Phi();
529 return c.Pt();
530}
531//_____________________________________________________________________________
532template <typename T> Float_t AliAnalysisTwoParticleResonanceFlowTask::PtSelector(Int_t tracktype, const T* track1, const T* track2, Float_t mass) const
533{
534 // plot m_inv spectra of like- and unlike sign pairs
535 Float_t pt = PairPt(track1, track2);
536 if (tracktype == 0) { // signal histograms
537 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
538 if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) {
539 fResonanceSignal[ptbin]->Fill(mass);
540 if(fQA) fPtSpectra[ptbin]->Fill(pt);
541 }
542 }
543 }
544 if (tracktype == 1) for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) if ((fPtBins[ptbin] <= pt) && (fPtBins[ptbin+1] > pt )) fResonanceBackground[ptbin]->Fill(mass);
545 return pt;
546}
547//_____________________________________________________________________________
548template <typename T> Bool_t AliAnalysisTwoParticleResonanceFlowTask::QualityCheck(T* track) const
549{
550 // check if track meets quality cuts
551 if(!track) return kFALSE;
552 return fPOICuts->IsSelected(track);
553}
554//_____________________________________________________________________________
8da4b0c9 555void AliAnalysisTwoParticleResonanceFlowTask::TrackQA(AliAODTrack* track, Int_t species, Bool_t allChargedParticles) const
556{
557 // fill qa histo's
558 if(allChargedParticles) { // track didn't pass identification yet
559 fNOPID->Fill(track->P(), track->GetTPCsignal());
560 if(track->Charge() > 0) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
561 if(track->Charge() < 0) {fEventStats->Fill(2); fPtN->Fill(track->Pt());}
562 }
563 else if(!allChargedParticles) { // track has been identified
564 if(species==fSpeciesA) {
565 fPtSpeciesA->Fill(track->Pt());
566 fPIDk->Fill(track->P(), track->GetTPCsignal());
567 }
568 else if(species==fSpeciesB) {
569 fPtSpeciesB->Fill(track->Pt());
570 fPIDp->Fill(track->P(), track->GetTPCsignal());
571 }
572 fPhi->Fill(track->Phi());
573 fEta->Fill(track->Eta());
574 }
575}
576//_____________________________________________________________________________
fbd0da9c 577template <typename T> void AliAnalysisTwoParticleResonanceFlowTask::SetNullCuts(T* event)
578{
579 // Set null cuts
fbd0da9c 580 fCutsRP->SetEvent(event, MCEvent());
9e437f94 581 if(event) fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
fbd0da9c 582 fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
583 fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
584 fNullCuts->SetEvent(event, MCEvent());
585}
586//_____________________________________________________________________________
587void AliAnalysisTwoParticleResonanceFlowTask::PrepareFlowEvent(Int_t iMulti)
588{
589 // Prepare flow events
fbd0da9c 590 fFlowEvent->ClearFast();
591 fFlowEvent->Fill(fCutsRP, fNullCuts);
592 fFlowEvent->SetReferenceMultiplicity(iMulti);
593 fFlowEvent->DefineDeadZone(0, 0, 0, 0);
594}
595//_____________________________________________________________________________
596void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethod(TObjArray* MixingCandidates)
597{
598 // extract v2 using the phi minus psi method
fbd0da9c 599 if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
fbd0da9c 600 return;
601 }
602 // retrieve data
603 Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
604 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
605 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
606 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
fbd0da9c 607 // if event plane stuff went ok, fill the histograms necessary for flow analysis with phi - psi method
fbd0da9c 608 AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
70adca5c 609 if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f, fatal error! ", fCentrality, fVertex));
610 else if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
fbd0da9c 611 Int_t nEvents = pool->GetCurrentNEvents();
fbd0da9c 612 for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
613 TObjArray* mixed_candidates = pool->GetEvent(iEvent);
614 if(!mixed_candidates) continue; // this should NEVER happen
615 Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
616 Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
fbd0da9c 617 TObjArray* SpeciesA = new TObjArray();
618 SpeciesA->SetOwner(kTRUE);
619 TObjArray* SpeciesB = new TObjArray();
620 SpeciesB->SetOwner(kTRUE);
621 TObjArray* SpeciesAFromBuffer = new TObjArray();
622 SpeciesAFromBuffer->SetOwner(kTRUE);
623 TObjArray* SpeciesBFromBuffer = new TObjArray();
624 SpeciesBFromBuffer->SetOwner(kTRUE);
625 for (Int_t iTracks = 0; iTracks < candidates; iTracks++) {
626 AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
627 if (!track) continue;
628 if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
629 else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
630 }
631 for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) {
632 AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
633 if (!track) continue;
634 if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
635 else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
636 }
637 PhiMinusPsiMethodWriteData(kTRUE, SpeciesA, SpeciesB, abcPsi2); // write signal information to histograms
638 PhiMinusPsiMethodWriteData(kFALSE, SpeciesA, SpeciesBFromBuffer, abcPsi2);
639 PhiMinusPsiMethodWriteData(kFALSE, SpeciesAFromBuffer, SpeciesB, abcPsi2);
640 } // end of mixed events loop
641 } // end of checking to see whether pool is filled correctly
642 pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
fbd0da9c 643}
644//_____________________________________________________________________________
645void AliAnalysisTwoParticleResonanceFlowTask::PhiMinusPsiMethodWriteData(Bool_t signal, TObjArray* SpeciesA, TObjArray* SpeciesB, Float_t* abcPsi2)
646{
647 // write data for phi minus psi method into correct histograms
fbd0da9c 648 for (Int_t i = 0; i < SpeciesA->GetEntriesFast(); i++) { // create pairs
649 for(Int_t j = 0; j < SpeciesB->GetEntriesFast(); j++) {
650 AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)(SpeciesA->At(i));
651 AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)(SpeciesB->At(j));
70adca5c 652 if(!(trackA&&trackB)) continue; // shouldn't happen
fbd0da9c 653 if(trackA->ID() == trackB->ID()) continue; // remove autocorrelations
654 if(fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
655 if(fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
656 Float_t minv = InvariantMass(trackA, trackB);
657 Float_t pt = PtSelector(2, trackA, trackB, minv);
658 TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
659 TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
660 TVector3 c = a + b;
661 Float_t phi = c.Phi();
662// if(phi < 0) phi += TMath::Pi();
663 Float_t dPhi_a = phi - abcPsi2[0];
664 Float_t dPhi_c = phi - abcPsi2[2];
665// if (dPhi_a < 0) dPhi_a += TMath::Pi();
666// if (dPhi_c < 0) dPhi_c += TMath::Pi();
667 // now all necessary numbers are here, so we can put them into histograms
668 for(Int_t k(0); k < fNdPhiBins; k++) { // dPhi bin loop
669 if(dPhi_a >= fdPhiBins[k] && dPhi_a < fdPhiBins[k+1]) // see if track is in desired dPhi_a range
670 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
671 if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
672 signal ? fPhiMinusPsiDataContainer[k][ptbin][0]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][0]->Fill(minv);
673 }
674 if(dPhi_c >= fdPhiBins[k] && dPhi_c < fdPhiBins[k+1]) // see if track is in desired dPhi_c range
675 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
676 if(pt <= fPtBins[ptbin] || pt > fPtBins[ptbin+1]) continue;
677 signal ? fPhiMinusPsiDataContainer[k][ptbin][1]->Fill(minv) : fPhiMinusPsiBackgroundContainer[k][ptbin][1]->Fill(minv);
678 }
679 }
680 } // end of loop on i- tracks
681 } // end of loop on i+ tracks
682}
683//_____________________________________________________________________________
684void AliAnalysisTwoParticleResonanceFlowTask::VZEROSubEventAnalysis()
685{
686 // vzero event plane analysis using three subevents
fbd0da9c 687 if(!AliAnalysisTaskVnV0::IsPsiComputed()) { // AliAnalysisTaskVnV0 must be added to analysis que before this task !!!
fbd0da9c 688 return;
689 }
690 // retrieve data
691 Float_t abcPsi2[] = {AliAnalysisTaskVnV0::GetPsi2V0A(), AliAnalysisTaskVnV0::GetPsi2TPC(), AliAnalysisTaskVnV0::GetPsi2V0C()};
692 for(Int_t i(0); i < 3; i++) if(abcPsi2[i] == 0) {
fbd0da9c 693 return;
694 }
695 // save info for resolution calculations
696 fSubEventDPhiv2->Fill(0.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[1]))); // vzeroa - tpc
697 fSubEventDPhiv2->Fill(1.5, TMath::Cos(2.*(abcPsi2[0]-abcPsi2[2]))); // vzeroa - vzeroc
698 fSubEventDPhiv2->Fill(2.5, TMath::Cos(2.*(abcPsi2[1]-abcPsi2[2]))); // tpc - vzeroc
699 Float_t minv[fMassBins+1];
700 Float_t _dphi[fMassBins][fNPtBins][2]; //minv, pt, v0a-c
701 Int_t occurence[fMassBins][fNPtBins]; //minv, pt
702 Float_t _inc = (fMaxMass - fMinMass) / (Float_t)fMassBins;
703 if(_inc==0) return; // avoid division by zero later on
704 for(Int_t mb(0); mb < fMassBins+1; mb++) minv[mb] = fMinMass + mb * _inc;
705 for(Int_t i(0); i < fMassBins; i++) for (Int_t j(0); j < fNPtBins; j++) for(Int_t k(0); k < 2; k ++) {
706 _dphi[i][j][k] = 0;
707 if(k==0) occurence[i][j] = 0;
708 }
709 for(Int_t iCand(0); iCand < fCandidates->GetEntriesFast(); iCand++) { // loop over unlike sign tracks
710 AliFlowTrackSimple *track = dynamic_cast<AliFlowTrackSimple*>(fCandidates->At(iCand));
711 if(!track) {
fbd0da9c 712 continue;
713 }
714 for(Int_t mb(0); mb < fMassBins; mb++) { // loop over massbands
715 if(track->Mass() <= minv[mb] || track->Mass() > minv[mb+1]) continue;
716 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) { // loop over pt bins
717 if(track->Pt() <= fPtBins[ptbin] || track->Pt() > fPtBins[ptbin+1]) continue;
718 _dphi[mb][ptbin][0]+=TMath::Cos(2.*(track->Phi() - abcPsi2[0]));
719 _dphi[mb][ptbin][1]+=TMath::Cos(2.*(track->Phi() - abcPsi2[2]));
720 occurence[mb][ptbin]+=1;
721 }
722 }
723 for(Int_t mb(0); mb < fMassBins; mb++) // write vn values to tprofiles
724 for(Int_t ptbin(0); ptbin < fNPtBins; ptbin++) {
725 if(occurence[mb][ptbin]==0) continue;
726 fV0Data[ptbin][0]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][0]/(Float_t)occurence[mb][ptbin]);
727 fV0Data[ptbin][1]->Fill(mb*_inc+fMinMass+0.5*_inc, _dphi[mb][ptbin][1]/(Float_t)occurence[mb][ptbin]);
728 }
729 }
730}
731//_____________________________________________________________________________
9e437f94 732void AliAnalysisTwoParticleResonanceFlowTask::DoAnalysisOnTheFly(AliFlowEventSimple* event)
fbd0da9c 733{
9e437f94 734 // initialize the task for on the fly analysis and call the user exec
735 if(fFlowEvent) delete fFlowEvent; // clear out the old flow event
736 fFlowEvent = (AliFlowEvent*)event; // cast the input event to its derived type
737 UserExec("fly"); // call the UserExec with flag 'on the fly'
738}
739//_____________________________________________________________________________
740void AliAnalysisTwoParticleResonanceFlowTask::DoAnalysisOnTheFly(TDirectoryFile* outputFile)
741{
742 // write the anlaysis to an output file
743 outputFile->Add(fOutputList);
744 outputFile->Write(outputFile->GetName(), TObject::kSingleKey);
745}
746//_____________________________________________________________________________
747void AliAnalysisTwoParticleResonanceFlowTask::DoAnalysisOnTheFly(TObjArray* MixingCandidates, TObjArray* SpeciesA, TObjArray* ocSpeciesA, TObjArray* SpeciesB, TObjArray* ocSpeciesB)
748{
749 // do the flow analysis on the fly.
750 Int_t iTracks = fFlowEvent->NumberOfTracks();
751 fCandidates->SetLast(-1); // clean out candidate array
752 Int_t charge(-1), tID(0); // we'll use this guy to store charge
753 if(fQA) fEventStats->Fill(0); // event counter
754 for (Int_t i = 0; i < iTracks; i++) { // track loop. iterator i is used as unique track id (necessary later on to avoid auto-correlations)
755 AliFlowTrackSimple* track = fFlowEvent->GetTrack(i);
471b00a8 756 if(!track) continue;
9e437f94 757 tID = track->GetID(); // store ID
758 (tID > 0) ? charge = 1 : charge = -1 ; // get the charge of the track
26f043a1 759 Double_t pt(track->Pt()), phi(track->Phi()), px(pt*TMath::Cos(phi)), py(pt*TMath::Sin(phi)), pz(track->Weight()), p(TMath::Sqrt(px*px+py*py+pz*pz)); // TODO ugly, but pz is stored as weight ...
9e437f94 760 if (charge == fChargeA && TMath::Abs(tID)==TMath::Abs(fSpeciesA)) { // store species a
26f043a1 761 SpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassA, i, fSpeciesA));
762 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassA, i, fSpeciesA));
9e437f94 763 }
764 if (charge == -1*fChargeA && TMath::Abs(tID)==TMath::Abs(fSpeciesA)) { // store opposite charge species a
26f043a1 765 ocSpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassA, i,fSpeciesA));
766 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassA, i, fSpeciesA));
9e437f94 767 }
768 if (charge == fChargeB && TMath::Abs(tID)==TMath::Abs(fSpeciesB)) { // store species b
26f043a1 769 SpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassB, i, fSpeciesB));
770 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassB, i, fSpeciesB));
9e437f94 771 }
772 if (charge == -1*fChargeB && TMath::Abs(tID)==TMath::Abs(fSpeciesB)) { // store opposite charge species b
26f043a1 773 ocSpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassB, i, fSpeciesB));
774 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), phi, p, px, py, pz, pt, charge, fMassB, i, fSpeciesB));
9e437f94 775 }
776 // at the end: convert the flow track to an 'actual' flow track with id and charge
777 track->SetID(i); // set the unique id
778 track->SetCharge(charge); // set charge
779 }
780}
781//_____________________________________________________________________________
782void AliAnalysisTwoParticleResonanceFlowTask::UserExec(Option_t * option)
783{
784 // UserExec: called for each event. Commented where necessary
fbd0da9c 785 TObjArray* MixingCandidates = 0x0; // init null pointer for event mixing
70adca5c 786 if(fEventMixing || fPhiMinusPsiMethod) {
fbd0da9c 787 MixingCandidates = new TObjArray();
788 MixingCandidates->SetOwner(kTRUE); // mixing candidates has ownership of objects in array
789 }
fbd0da9c 790 TObjArray* SpeciesA = new TObjArray(); // create arrays for the helper tracks
791 SpeciesA->SetOwner(kTRUE);
792 TObjArray* ocSpeciesA = new TObjArray(); // opposite charge particles
793 ocSpeciesA->SetOwner(kTRUE);
794 TObjArray* SpeciesB = new TObjArray();
795 SpeciesB->SetOwner(kTRUE);
796 TObjArray* ocSpeciesB = new TObjArray();
797 ocSpeciesB->SetOwner(kTRUE);
9e437f94 798 // check the options
799 char chopt[128];
800 strlcpy(chopt,option,128);
801 char *onTheFly;
802 onTheFly = strstr(chopt,"fly");
803 if(onTheFly) { // do the on the fly analysis
804 printf(" > we're ready to fly ... ! \n");
70adca5c 805 MixingCandidates = new TObjArray();
806 MixingCandidates->SetOwner(kTRUE);
9e437f94 807 DoAnalysisOnTheFly(MixingCandidates, SpeciesA, ocSpeciesA, SpeciesB, ocSpeciesB);
fbd0da9c 808 }
9e437f94 809 if(!onTheFly) { // do analysis in the regular way (with the manager, etc)
810 if (!fPIDResponse) { // kill the event if there isn't a pid response
811 return;
fbd0da9c 812 }
9e437f94 813
814 fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
815 if (fAOD) {
f5cd76c4 816 if (fCutsEvent && !fCutsEvent->PassesCuts(fAOD, 0x0)) return; // check for event cuts
9e437f94 817 InitializeBayesianPID(fAOD); // init event objects
818 SetNullCuts(fAOD);
819 Int_t iTracks = fAOD->GetNumberOfTracks();
820 PrepareFlowEvent(iTracks); // does number of tracks correspond to the set filterbit ??? FIXME !!!
821 fCandidates->SetLast(-1);
822 if(fIsMC) IsMC(); // launch mc stuff FIXME
823 if(fQA) fEventStats->Fill(0);
824 for (Int_t i = 0; i < iTracks; i++) { // select analysis candidates
825 AliAODTrack* track = fAOD->GetTrack(i);
826 if (!QualityCheck(track)) continue; // reject poor quality tracks
827 if (track->Charge() == fChargeA && AcceptTrack(track, fSpeciesA)) { // store species a
828 SpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
829 track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
830 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
831 track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA,
832 track->GetID(), fSpeciesA));
833 }
834 if (track->Charge() == -1*fChargeA && AcceptTrack(track, fSpeciesA)) { // store opposite charge species a
835 ocSpeciesA->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
836 track->Pt(), track->Charge(), fMassA, track->GetID(), fSpeciesA));
837 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
838 track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassA,
839 track->GetID(), fSpeciesA));
840 }
841 if (track->Charge() == fChargeB && AcceptTrack(track, fSpeciesB)) { // store species b
842 SpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
843 track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
844 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
845 track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB,
846 track->GetID(), fSpeciesB));
847 }
848 if (track->Charge() == -1*fChargeB && AcceptTrack(track, fSpeciesB)) { // store opposite charge species b
849 ocSpeciesB->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(), track->Py(), track->Pz(),
850 track->Pt(), track->Charge(), fMassB, track->GetID(), fSpeciesB));
851 if(fEventMixing) MixingCandidates->Add(new AliResonanceFlowHelperTrack(track->Eta(), track->Phi(), track->P(), track->Px(),
852 track->Py(), track->Pz(), track->Pt(), track->Charge(), fMassB,
853 track->GetID(), fSpeciesB));
854
855 }
856 }
857 } // end the loop for event manager's aod's
858 } // end of data loop
859 // do the phi minus psi method if that's specified FIXME currenlty only supports event mixing
860 if (fPhiMinusPsiMethod) { // for the phi minus psi method, no need for flow package etc
861 PhiMinusPsiMethod(MixingCandidates); // call the method
862 PostData(1, fOutputList); // save the output data
863 return; // return to retrieve next event
864 }
865 // start the resonance reconstruction, this fills the fCandidates array with flow tracks
866 ResonanceSignal(SpeciesA, SpeciesB);
867 if(fV0) VZEROSubEventAnalysis();
868 for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
869 AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
870 if (!cand) continue;
871 for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
872 for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
873 AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
874 if (!iRP) continue;
875 if (!iRP->InRPSelection()) continue;
876 if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
877 iRP->SetForRPSelection(kFALSE);
878 fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
fbd0da9c 879 }
fbd0da9c 880 }
fbd0da9c 881 }
9e437f94 882 cand->SetForPOISelection(kTRUE);
883 fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
884 }
885 if(!fEventMixing) { // do the combinatorial background
886 ResonanceBackground(ocSpeciesA, SpeciesB); // mix opposite charge species A with species B
887 if(fSpeciesA!=fSpeciesB) ResonanceBackground(SpeciesA, ocSpeciesB); // mix species A with opposite charge species B
888 }
889 else ReconstructionWithEventMixing(MixingCandidates);
890 if(!onTheFly) {
891 PostData(1, fOutputList);
892 PostData(2, fFlowEvent);
fbd0da9c 893 }
9e437f94 894 delete SpeciesA;
895 delete SpeciesB;
896 delete ocSpeciesA;
897 delete ocSpeciesB; // clean heap memory
fbd0da9c 898}
899//_____________________________________________________________________________
900void AliAnalysisTwoParticleResonanceFlowTask::ResonanceSignal(TObjArray* SpeciesA, TObjArray* SpeciesB) const
901{
902 // fill signal histograms
9e437f94 903 Int_t spA(SpeciesA->GetEntries()), spB(SpeciesB->GetEntries());
904 for (Int_t i(0); i < spA; i++) { //track loop over species A
905 for (Int_t j(0); j < spB; j++) { //track loop over species B
fbd0da9c 906 AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
907 AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
70adca5c 908 if(!(trackA&&trackB)) continue; // shouldn't happen
fbd0da9c 909 if(trackA->ID()==trackB->ID()) continue; // beware of autocorrelations
910 if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
911 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
912 Float_t mass = InvariantMass(trackA, trackB);
913 Float_t pt = PtSelector(0, trackA, trackB, mass);
914 TVector3 a(trackA->Px(), trackA->Py(), trackA->Pz());
915 TVector3 b(trackB->Px(), trackB->Py(), trackB->Pz());
916 TVector3 c = a + b;
917 Float_t phi = c.Phi();
918 Float_t eta = c.Eta();
919 Int_t nIDs[2];
920 nIDs[0] = trackA->ID();
921 nIDs[1] = trackB->ID();
922 MakeTrack(mass, pt, phi, eta, 2, nIDs);
923 }
924 }
925}
926//_____________________________________________________________________________
927void AliAnalysisTwoParticleResonanceFlowTask::ResonanceBackground(TObjArray* SpeciesA, TObjArray* SpeciesB, Bool_t checkAutoCorrelations) const
928{
929 // fill background histograms
930 for (Int_t i(0); i < SpeciesA->GetEntriesFast(); i++) { //track loop over species A
931 for (Int_t j(0); j < SpeciesB->GetEntriesFast(); j++) { //track loop over species B
932 AliResonanceFlowHelperTrack* trackA = (AliResonanceFlowHelperTrack*)SpeciesA->At(i);
933 AliResonanceFlowHelperTrack* trackB = (AliResonanceFlowHelperTrack*)SpeciesB->At(j);
70adca5c 934 if(!(trackA&&trackB)) continue; // shouldn't happen
fbd0da9c 935 if((trackA->ID() == trackB->ID()) && checkAutoCorrelations) continue; // remove autocorrelations
936 if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(trackA, trackB))) continue;
937 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(trackA, trackB))) continue;
938 PtSelector(1, trackA, trackB, InvariantMass(trackA, trackB));
939 }
940 }
941}
942//_____________________________________________________________________________
943void AliAnalysisTwoParticleResonanceFlowTask::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
944{
945 // perform reconstruction with event mixing
fbd0da9c 946 AliEventPool* pool = fPoolManager->GetEventPool(fCentrality, fVertex);
70adca5c 947 if(!pool) AliFatal(Form("No pool found for centrality = %f, zVtx = %f, fatal error!", fCentrality, fVertex));
948 else if(pool->IsReady() || pool->NTracksInPool() > fMixingParameters[1] / 10 || pool->GetCurrentNEvents() >= fMixingParameters[2]) {
fbd0da9c 949 Int_t nEvents = pool->GetCurrentNEvents();
fbd0da9c 950 for (Int_t iEvent(0); iEvent < nEvents; iEvent++) {
951 TObjArray* mixed_candidates = pool->GetEvent(iEvent);
952 if(!mixed_candidates) continue; // this should NEVER happen
953 Int_t bufferTracks = mixed_candidates->GetEntriesFast(); // buffered candidates
954 Int_t candidates = MixingCandidates->GetEntriesFast(); // mixing candidates
fbd0da9c 955 TObjArray* SpeciesA = new TObjArray();
956 SpeciesA->SetOwner(kTRUE);
957 TObjArray* SpeciesB = new TObjArray();
958 SpeciesB->SetOwner(kTRUE);
959 TObjArray* SpeciesAFromBuffer = new TObjArray();
960 SpeciesAFromBuffer->SetOwner(kTRUE);
961 TObjArray* SpeciesBFromBuffer = new TObjArray();
962 SpeciesBFromBuffer->SetOwner(kTRUE);
963 for (Int_t iTracks = 0; iTracks < candidates; iTracks++) {
964 AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)MixingCandidates->At(iTracks);
965 if (!track) continue;
966 if (track->Charge() == fChargeA && track->Species() == fSpeciesA) SpeciesA->Add(track);
967 else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesB->Add(track);
968 }
969 for (Int_t iTracks = 0; iTracks < bufferTracks; iTracks++) {
970 AliResonanceFlowHelperTrack* track = (AliResonanceFlowHelperTrack*)mixed_candidates->At(iTracks);
971 if (!track) continue;
972 if (track->Charge() == fChargeA && track->Species() == fSpeciesA ) SpeciesAFromBuffer->Add(track);
973 else if (track->Charge() == fChargeB && track->Species() == fSpeciesB) SpeciesBFromBuffer->Add(track);
974 }
975 ResonanceBackground(SpeciesA, SpeciesBFromBuffer, kFALSE);
976 ResonanceBackground(SpeciesB, SpeciesAFromBuffer, kFALSE);
977 } // end of mixed events loop
978 } // end of checking to see whether pool is filled correctly
979 pool->UpdatePool(MixingCandidates); // update pool with current mixing candidates (for next event)
fbd0da9c 980}
981//_____________________________________________________________________________
982void AliAnalysisTwoParticleResonanceFlowTask::Terminate(Option_t *)
983{
984 // terminate
fbd0da9c 985}
986//______________________________________________________________________________
987void AliAnalysisTwoParticleResonanceFlowTask::MakeTrack(Float_t mass, Float_t pt, Float_t phi, Float_t eta, Int_t nDau, Int_t iID[]) const
988{
989 // Construct Flow Candidate Track from two selected candidates
fbd0da9c 990 Bool_t overwrite = kTRUE;
991 AliFlowCandidateTrack *sTrack = static_cast<AliFlowCandidateTrack*>(fCandidates->At(fCandidates->GetLast() + 1));
992 if (!sTrack) {
993 sTrack = new AliFlowCandidateTrack(); //deleted by fCandidates
994 overwrite = kFALSE;
995 }
996 else sTrack->ClearMe();
997 sTrack->SetMass(mass);
998 sTrack->SetPt(pt);
999 sTrack->SetPhi(phi);
1000 sTrack->SetEta(eta);
1001 for (Int_t iDau = 0; iDau != nDau; ++iDau) sTrack->AddDaughter(iID[iDau]);
1002 sTrack->SetForPOISelection(kTRUE);
1003 sTrack->SetForRPSelection(kFALSE);
1004 if (overwrite) fCandidates->SetLast(fCandidates->GetLast() + 1);
1005 else fCandidates->AddLast(sTrack);
1006 return;
1007}
1008//_____________________________________________________________________________
1009void AliAnalysisTwoParticleResonanceFlowTask::IsMC()
1010{
1011 ForceExit(1,"No MC mode available at this stage of developement ... ");
1012}
1013//_____________________________________________________________________________