]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/FLOW/Tasks/AliAnalysisTaskVnV0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
some coverity fixes
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskVnV0.cxx
CommitLineData
afa8df58 1#include "AliAnalysisTaskVnV0.h"
2
3// ROOT includes
4#include <TMath.h>
5
6// AliRoot includes
7#include "AliInputEventHandler.h"
8#include "AliAODEvent.h"
9#include "AliAODVertex.h"
10#include "AliAODTrack.h"
11#include "AliCentrality.h"
12#include "AliVHeader.h"
13#include "AliAODVZERO.h"
14#include "TFile.h"
15#include "AliOADBContainer.h"
16#include "TH2F.h"
17#include "TF1.h"
587d006a 18#include "AliGenHijingEventHeader.h"
19#include "AliMCEvent.h"
20#include "AliAODMCHeader.h"
21#include "AliAODMCParticle.h"
243fbce7 22#include "TChain.h"
a640d925 23#include "AliESDtrackCuts.h"
24#include "AliESDVertex.h"
a31d07c8 25#include "AliEventplane.h"
26#include "TProfile2D.h"
afa8df58 27
28// STL includes
29//#include <iostream>
30//using namespace std;
31
32ClassImp(AliAnalysisTaskVnV0)
6fa5bb6e 33Bool_t AliAnalysisTaskVnV0::fgIsPsiComputed = kFALSE;
34Float_t AliAnalysisTaskVnV0::fgPsi2v0a=999.;
35Float_t AliAnalysisTaskVnV0::fgPsi2v0c=999.;
36Float_t AliAnalysisTaskVnV0::fgPsi2tpc=999.;
37Float_t AliAnalysisTaskVnV0::fgPsi3v0a=999.;
38Float_t AliAnalysisTaskVnV0::fgPsi3v0c=999.;
39Float_t AliAnalysisTaskVnV0::fgPsi3tpc=999.;
a640d925 40Float_t AliAnalysisTaskVnV0::fgPsi2v0aMC=999.;
41Float_t AliAnalysisTaskVnV0::fgPsi2v0cMC=999.;
42Float_t AliAnalysisTaskVnV0::fgPsi2tpcMC=999.;
43Float_t AliAnalysisTaskVnV0::fgPsi3v0aMC=999.;
44Float_t AliAnalysisTaskVnV0::fgPsi3v0cMC=999.;
45Float_t AliAnalysisTaskVnV0::fgPsi3tpcMC=999.;
afa8df58 46
afa8df58 47//_____________________________________________________________________________
48AliAnalysisTaskVnV0::AliAnalysisTaskVnV0():
49 AliAnalysisTaskSE(),
587d006a 50 fVtxCut(10.0), // cut on |vertex| < fVtxCut
afa8df58 51 fEtaCut(0.8), // cut on |eta| < fEtaCut
52 fMinPt(0.15), // cut on pt > fMinPt
a640d925 53 fMinDistV0(0),
54 fMaxDistV0(100),
243fbce7 55 fV2(kTRUE),
56 fV3(kTRUE),
57 fIsMC(kFALSE),
58 fQAsw(kFALSE),
a31d07c8 59 fIsAfter2011(kFALSE),
afa8df58 60 fRun(-1),
a640d925 61 fNcluster(70),
afa8df58 62 fList(new TList()),
63 fList2(new TList()),
587d006a 64 fList3(new TList()),
65 fList4(new TList()),
66 fMultV0(NULL),
afa8df58 67 fV0Cpol(100),
68 fV0Apol(100),
587d006a 69 fHResTPCv0A2(NULL),
70 fHResTPCv0C2(NULL),
71 fHResv0Cv0A2(NULL),
72 fHResTPCv0A3(NULL),
73 fHResTPCv0C3(NULL),
74 fHResv0Cv0A3(NULL),
75 fPhiRPv0A(NULL),
76 fPhiRPv0C(NULL),
77 fPhiRPv0Av3(NULL),
78 fPhiRPv0Cv3(NULL),
587d006a 79 fQA(NULL),
80 fQA2(NULL),
81 fQAv3(NULL),
82 fQA2v3(NULL),
afa8df58 83 fPID(new AliFlowBayesianPID()),
587d006a 84 fTree(NULL),
85 fCentrality(-1),
afa8df58 86 evPlAngV0ACor2(0),
87 evPlAngV0CCor2(0),
88 evPlAng2(0),
89 evPlAngV0ACor3(0),
90 evPlAngV0CCor3(0),
91 evPlAng3(0),
587d006a 92 fContAllChargesV0A(NULL),
93 fContAllChargesV0C(NULL),
94 fContAllChargesV0Av3(NULL),
95 fContAllChargesV0Cv3(NULL),
96 fContAllChargesMC(NULL),
08371bdc 97 fHResMA2(NULL),
98 fHResMC2(NULL),
99 fHResAC2(NULL),
100 fHResMA3(NULL),
101 fHResMC3(NULL),
102 fHResAC3(NULL),
103 fContAllChargesMCA(NULL),
104 fContAllChargesMCC(NULL),
105 fContAllChargesMCAv3(NULL),
0f25ad32 106 fContAllChargesMCCv3(NULL),
107 fFillDCA(kFALSE),
108 fContQApid(NULL),
a640d925 109 fModulationDEDx(kFALSE),
a31d07c8 110 fZvtx(0.),
111 fNK0s(0),
112 fNpiPos(0),
113 fNpiNeg(0),
114 fHKsPhi(NULL),
115 fHKsPhiEP(NULL),
116 fHK0sMass(NULL),
117 fHK0sMass2(NULL),
118 fHK0vsLambda(NULL),
119 fHctauPtEP(NULL),
120 fHctauAt1EP(NULL),
a640d925 121 fCutsDaughter(NULL)
afa8df58 122{
afa8df58 123 // Default constructor (should not be used)
124 fList->SetName("resultsV2");
125 fList2->SetName("resultsV3");
587d006a 126 fList3->SetName("resultsMC");
127 fList4->SetName("QA");
afa8df58 128
243fbce7 129 fList->SetOwner(kTRUE);
130 fList2->SetOwner(kTRUE);
131 fList3->SetOwner(kTRUE);
132 fList4->SetOwner(kTRUE);
133
afa8df58 134 fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
455f6523 135
136 for(Int_t i=0;i < 1000;i++){
137 fPhiK0s[i] = 0.0;
138 fPtK0s[i] = 0.0;
139 fIPiPos[i] = 0;
140 fIPiNeg[i] = 0;
141 }
afa8df58 142}
143
144//______________________________________________________________________________
145AliAnalysisTaskVnV0::AliAnalysisTaskVnV0(const char *name):
146 AliAnalysisTaskSE(name),
587d006a 147 fVtxCut(10.0), // cut on |vertex| < fVtxCut
afa8df58 148 fEtaCut(0.8), // cut on |eta| < fEtaCut
149 fMinPt(0.15), // cut on pt > fMinPt
a640d925 150 fMinDistV0(0),
151 fMaxDistV0(100),
243fbce7 152 fV2(kTRUE),
153 fV3(kTRUE),
154 fIsMC(kFALSE),
155 fQAsw(kFALSE),
a31d07c8 156 fIsAfter2011(kFALSE),
afa8df58 157 fRun(-1),
a640d925 158 fNcluster(70),
afa8df58 159 fList(new TList()),
160 fList2(new TList()),
587d006a 161 fList3(new TList()),
162 fList4(new TList()),
163 fMultV0(NULL),
afa8df58 164 fV0Cpol(100),
165 fV0Apol(100),
587d006a 166 fHResTPCv0A2(NULL),
167 fHResTPCv0C2(NULL),
168 fHResv0Cv0A2(NULL),
169 fHResTPCv0A3(NULL),
170 fHResTPCv0C3(NULL),
171 fHResv0Cv0A3(NULL),
172 fPhiRPv0A(NULL),
173 fPhiRPv0C(NULL),
174 fPhiRPv0Av3(NULL),
175 fPhiRPv0Cv3(NULL),
587d006a 176 fQA(NULL),
177 fQA2(NULL),
178 fQAv3(NULL),
179 fQA2v3(NULL),
afa8df58 180 fPID(new AliFlowBayesianPID()),
587d006a 181 fTree(NULL),
182 fCentrality(-1),
afa8df58 183 evPlAngV0ACor2(0),
184 evPlAngV0CCor2(0),
185 evPlAng2(0),
186 evPlAngV0ACor3(0),
187 evPlAngV0CCor3(0),
188 evPlAng3(0),
587d006a 189 fContAllChargesV0A(NULL),
190 fContAllChargesV0C(NULL),
191 fContAllChargesV0Av3(NULL),
192 fContAllChargesV0Cv3(NULL),
193 fContAllChargesMC(NULL),
08371bdc 194 fHResMA2(NULL),
195 fHResMC2(NULL),
196 fHResAC2(NULL),
197 fHResMA3(NULL),
198 fHResMC3(NULL),
199 fHResAC3(NULL),
200 fContAllChargesMCA(NULL),
201 fContAllChargesMCC(NULL),
202 fContAllChargesMCAv3(NULL),
0f25ad32 203 fContAllChargesMCCv3(NULL),
204 fFillDCA(kFALSE),
205 fContQApid(NULL),
a640d925 206 fModulationDEDx(kFALSE),
a31d07c8 207 fZvtx(0.),
208 fNK0s(0),
209 fNpiPos(0),
210 fNpiNeg(0),
211 fHKsPhi(NULL),
212 fHKsPhiEP(NULL),
213 fHK0sMass(NULL),
214 fHK0sMass2(NULL),
215 fHK0vsLambda(NULL),
216 fHctauPtEP(NULL),
217 fHctauAt1EP(NULL),
a640d925 218 fCutsDaughter(NULL)
afa8df58 219{
220
221 DefineOutput(1, TList::Class());
222 DefineOutput(2, TList::Class());
587d006a 223 DefineOutput(3, TList::Class());
224 DefineOutput(4, TList::Class());
afa8df58 225
226 // Output slot #1 writes into a TTree
227 fList->SetName("resultsV2");
228 fList2->SetName("resultsV3");
587d006a 229 fList3->SetName("resultsMC");
230 fList4->SetName("QA");
afa8df58 231
243fbce7 232 fList->SetOwner(kTRUE);
233 fList2->SetOwner(kTRUE);
234 fList3->SetOwner(kTRUE);
235 fList4->SetOwner(kTRUE);
236
afa8df58 237 fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
455f6523 238
239 for(Int_t i=0;i < 1000;i++){
240 fPhiK0s[i] = 0.0;
241 fPtK0s[i] = 0.0;
242 fIPiPos[i] = 0;
243 fIPiNeg[i] = 0;
244 }
afa8df58 245}
afa8df58 246//_____________________________________________________________________________
247AliAnalysisTaskVnV0::~AliAnalysisTaskVnV0()
248{
249
250}
251
252//______________________________________________________________________________
253void AliAnalysisTaskVnV0::UserCreateOutputObjects()
254{
255
587d006a 256 if(fIsMC) fPID->SetMC(kTRUE);
257
258
afa8df58 259 // Tree for EP debug (comment the adding to v2 list id not needed)
260 fTree = new TTree("tree","tree");
261 fTree->Branch("evPlAngV0ACor2",&evPlAngV0ACor2,"evPlAngV0ACor2/F");
262 fTree->Branch("evPlAngV0CCor2",&evPlAngV0CCor2,"evPlAngV0CCor2/F");
263 fTree->Branch("evPlAng2",&evPlAng2,"evPlAng2/F");
264 fTree->Branch("fCentrality",&fCentrality,"fCentrality/F");
265 fTree->Branch("evPlAngV0ACor3",&evPlAngV0ACor3,"evPlAngV0ACor3/F");
266 fTree->Branch("evPlAngV0CCor3",&evPlAngV0CCor3,"evPlAngV0CCor3/F");
267 fTree->Branch("evPlAng3",&evPlAng3,"evPlAng3/F");
268
269
270 // Container analyses (different steps mean different species)
271 const Int_t nPtBinsTOF = 45;
272 Double_t binsPtTOF[nPtBinsTOF+1] = {0., 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.25, 2.5, 2.75,3.0,3.25,3.5,3.75,4.0,4.5,5,5.5,6,6.5,7,8,9,10,12,15,20};
a640d925 273 const Int_t nCentrTOF = nCentrBin;
243fbce7 274 const Int_t nPsiTOF = 10;
587d006a 275 const Int_t nChargeBinsTOFres = 2;
a640d925 276 const Int_t nCentrTOFres = nCentrBin;
587d006a 277 const Int_t nProbTOFres = 4;
278 const Int_t nPsiTOFres = 10;
243fbce7 279 const Int_t nMaskPID = 3;
587d006a 280
0f25ad32 281 Int_t nDCABin = 1; // put to 1 not to store this info
282 if(fFillDCA) nDCABin = 3;
283 if(fIsMC && nDCABin>1) nDCABin = 6;
284 /*
285 0 = DCAxy < 2.4 && all (or Physical primary if MC)
286 1 = DCAxy > 2.4 && all (or Physical primary if MC)
287 2 = DCAxy < 2.4 && not Physical Primary for MC
288 3 = DCAxy > 2.4 && not Physical Primary for MC
289 */
290
291 Int_t binsTOF[6] = {nCentrTOFres,nChargeBinsTOFres,nProbTOFres,nPsiTOFres,nMaskPID,nDCABin};
587d006a 292 Int_t binsTOFmc[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,2};
08371bdc 293 Int_t binsTOFmcPureMC[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,1};
afa8df58 294
295 // v2 container
0f25ad32 296 fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
a640d925 297 fContAllChargesV0A->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
587d006a 298 fContAllChargesV0A->SetVarRange(1,-1.5,1.5); // charge
299 fContAllChargesV0A->SetVarRange(2,0.6,1.0001);// prob
300 fContAllChargesV0A->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
301 fContAllChargesV0A->SetVarRange(4,-0.5,2.5); // pid mask
0f25ad32 302 fContAllChargesV0A->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
587d006a 303 fContAllChargesV0A->SetVarName(0,"centrality");
304 fContAllChargesV0A->SetVarName(1,"charge");
305 fContAllChargesV0A->SetVarName(2,"prob");
306 fContAllChargesV0A->SetVarName(3,"#Psi");
307 fContAllChargesV0A->SetVarName(4,"PIDmask");
0f25ad32 308 fContAllChargesV0A->SetVarName(5,"DCAbin");
243fbce7 309 if(fV2) fContAllChargesV0A->AddSpecies("all",nPtBinsTOF,binsPtTOF);
310 if(fV2) fContAllChargesV0A->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
311 if(fV2) fContAllChargesV0A->AddSpecies("k",nPtBinsTOF,binsPtTOF);
312 if(fV2) fContAllChargesV0A->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
313 if(fV2) fContAllChargesV0A->AddSpecies("e",nPtBinsTOF,binsPtTOF);
314 if(fV2) fContAllChargesV0A->AddSpecies("d",nPtBinsTOF,binsPtTOF);
315 if(fV2) fContAllChargesV0A->AddSpecies("t",nPtBinsTOF,binsPtTOF);
316 if(fV2) fContAllChargesV0A->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
317 if(fV2) fContAllChargesV0A->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
a640d925 318 if(fV2) fContAllChargesV0A->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
319 if(fV2) fContAllChargesV0A->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
320 if(fV2) fContAllChargesV0A->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
a31d07c8 321 if(fV2) fContAllChargesV0A->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
587d006a 322
0f25ad32 323 fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
a640d925 324 fContAllChargesV0C->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
587d006a 325 fContAllChargesV0C->SetVarRange(1,-1.5,1.5); // charge
326 fContAllChargesV0C->SetVarRange(2,0.6,1.0001);// prob
327 fContAllChargesV0C->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
328 fContAllChargesV0C->SetVarRange(4,-0.5,2.5); // pid mask
0f25ad32 329 fContAllChargesV0C->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
587d006a 330 fContAllChargesV0C->SetVarName(0,"centrality");
331 fContAllChargesV0C->SetVarName(1,"charge");
332 fContAllChargesV0C->SetVarName(2,"prob");
333 fContAllChargesV0C->SetVarName(3,"#Psi");
334 fContAllChargesV0C->SetVarName(4,"PIDmask");
0f25ad32 335 fContAllChargesV0C->SetVarName(5,"DCAbin");
243fbce7 336 if(fV2) fContAllChargesV0C->AddSpecies("all",nPtBinsTOF,binsPtTOF);
337 if(fV2) fContAllChargesV0C->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
338 if(fV2) fContAllChargesV0C->AddSpecies("k",nPtBinsTOF,binsPtTOF);
339 if(fV2) fContAllChargesV0C->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
340 if(fV2) fContAllChargesV0C->AddSpecies("e",nPtBinsTOF,binsPtTOF);
341 if(fV2) fContAllChargesV0C->AddSpecies("d",nPtBinsTOF,binsPtTOF);
342 if(fV2) fContAllChargesV0C->AddSpecies("t",nPtBinsTOF,binsPtTOF);
343 if(fV2) fContAllChargesV0C->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
344 if(fV2) fContAllChargesV0C->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
a640d925 345 if(fV2) fContAllChargesV0C->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
346 if(fV2) fContAllChargesV0C->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
347 if(fV2) fContAllChargesV0C->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
a31d07c8 348 if(fV2) fContAllChargesV0C->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
afa8df58 349
350 fList->Add(fContAllChargesV0A);
351 fList->Add(fContAllChargesV0C);
352
a31d07c8 353 fHctauPtEP = new TProfile2D("hctauPtEP","K^{0}_{s} decay length;p_{T} (GeV/#it{c});#Delta#phi (rad)",40,0,5,10,-TMath::Pi(),TMath::Pi());
354 fHctauAt1EP = new TH2F("hctauAt1EP","K^{0}_{s} decay length at 1 GeV/#it{c};c#tau (cm);#Delta#phi (rad)",50,0,50,10,-TMath::Pi(),TMath::Pi());
355 // added at the end
356
243fbce7 357 if(fIsMC && fV2){
587d006a 358 fContAllChargesMC = new AliFlowVZEROResults("v2mc",5,binsTOFmc);
a640d925 359 fContAllChargesMC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
587d006a 360 fContAllChargesMC->SetVarRange(1,-1.5,1.5); // charge
361 fContAllChargesMC->SetVarRange(2,0.6,1.0001);// prob
362 fContAllChargesMC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
363 fContAllChargesMC->SetVarRange(4,-0.5,1.5); // pid mask
364 fContAllChargesMC->SetVarName(0,"centrality");
365 fContAllChargesMC->SetVarName(1,"charge");
366 fContAllChargesMC->SetVarName(2,"prob");
367 fContAllChargesMC->SetVarName(3,"#Psi");
368 fContAllChargesMC->SetVarName(4,"PIDmask");
369 fContAllChargesMC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
370 fContAllChargesMC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
371 fContAllChargesMC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
372 fContAllChargesMC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
373 fContAllChargesMC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
374 fContAllChargesMC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
375 fList3->Add(fContAllChargesMC);
08371bdc 376
377 fContAllChargesMCA = new AliFlowVZEROResults("v2mcA",5,binsTOFmcPureMC);
a640d925 378 fContAllChargesMCA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
08371bdc 379 fContAllChargesMCA->SetVarRange(1,-1.5,1.5); // charge
380 fContAllChargesMCA->SetVarRange(2,0.6,1.0001);// prob
381 fContAllChargesMCA->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
382 fContAllChargesMCA->SetVarRange(4,-0.5,1.5); // pid mask
383 fContAllChargesMCA->SetVarName(0,"centrality");
384 fContAllChargesMCA->SetVarName(1,"charge");
385 fContAllChargesMCA->SetVarName(2,"prob");
386 fContAllChargesMCA->SetVarName(3,"#Psi");
387 fContAllChargesMCA->SetVarName(4,"PIDmask");
388 fContAllChargesMCA->AddSpecies("all",nPtBinsTOF,binsPtTOF);
389 fContAllChargesMCA->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
390 fContAllChargesMCA->AddSpecies("k",nPtBinsTOF,binsPtTOF);
391 fContAllChargesMCA->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
392 fContAllChargesMCA->AddSpecies("e",nPtBinsTOF,binsPtTOF);
393 fContAllChargesMCA->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
394 fList3->Add(fContAllChargesMCA);
395
396 fContAllChargesMCC = new AliFlowVZEROResults("v2mcC",5,binsTOFmcPureMC);
a640d925 397 fContAllChargesMCC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
08371bdc 398 fContAllChargesMCC->SetVarRange(1,-1.5,1.5); // charge
399 fContAllChargesMCC->SetVarRange(2,0.6,1.0001);// prob
400 fContAllChargesMCC->SetVarRange(3,-TMath::Pi()/2,TMath::Pi()/2); // Psi
401 fContAllChargesMCC->SetVarRange(4,-0.5,1.5); // pid mask
402 fContAllChargesMCC->SetVarName(0,"centrality");
403 fContAllChargesMCC->SetVarName(1,"charge");
404 fContAllChargesMCC->SetVarName(2,"prob");
405 fContAllChargesMCC->SetVarName(3,"#Psi");
406 fContAllChargesMCC->SetVarName(4,"PIDmask");
407 fContAllChargesMCC->AddSpecies("all",nPtBinsTOF,binsPtTOF);
408 fContAllChargesMCC->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
409 fContAllChargesMCC->AddSpecies("k",nPtBinsTOF,binsPtTOF);
410 fContAllChargesMCC->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
411 fContAllChargesMCC->AddSpecies("e",nPtBinsTOF,binsPtTOF);
412 fContAllChargesMCC->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
413 fList3->Add(fContAllChargesMCC);
587d006a 414 }
415
afa8df58 416 // v3 container
0f25ad32 417 fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
a640d925 418 fContAllChargesV0Av3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
587d006a 419 fContAllChargesV0Av3->SetVarRange(1,-1.5,1.5); // charge
420 fContAllChargesV0Av3->SetVarRange(2,0.6,1.0001);// prob
421 fContAllChargesV0Av3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
0f25ad32 422 fContAllChargesV0Av3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
587d006a 423 fContAllChargesV0Av3->SetVarRange(4,-0.5,2.5); // pid mask
424 fContAllChargesV0Av3->SetVarName(0,"centrality");
425 fContAllChargesV0Av3->SetVarName(1,"charge");
426 fContAllChargesV0Av3->SetVarName(2,"prob");
427 fContAllChargesV0Av3->SetVarName(3,"#Psi");
428 fContAllChargesV0Av3->SetVarName(4,"PIDmask");
0f25ad32 429 fContAllChargesV0Av3->SetVarName(5,"DCAbin");
243fbce7 430 if(fV3) fContAllChargesV0Av3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
431 if(fV3) fContAllChargesV0Av3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
432 if(fV3) fContAllChargesV0Av3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
433 if(fV3) fContAllChargesV0Av3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
434 if(fV3) fContAllChargesV0Av3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
435 if(fV3) fContAllChargesV0Av3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
436 if(fV3) fContAllChargesV0Av3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
437 if(fV3) fContAllChargesV0Av3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
438 if(fV3) fContAllChargesV0Av3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
a640d925 439 if(fV3) fContAllChargesV0Av3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
440 if(fV3) fContAllChargesV0Av3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
441 if(fV3) fContAllChargesV0Av3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
a31d07c8 442 if(fV3) fContAllChargesV0Av3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
587d006a 443
0f25ad32 444 fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
a640d925 445 fContAllChargesV0Cv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
587d006a 446 fContAllChargesV0Cv3->SetVarRange(1,-1.5,1.5); // charge
447 fContAllChargesV0Cv3->SetVarRange(2,0.6,1.0001);// prob
448 fContAllChargesV0Cv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
449 fContAllChargesV0Cv3->SetVarRange(4,-0.5,2.5); // pid mask
0f25ad32 450 fContAllChargesV0Cv3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
587d006a 451 fContAllChargesV0Cv3->SetVarName(0,"centrality");
452 fContAllChargesV0Cv3->SetVarName(1,"charge");
453 fContAllChargesV0Cv3->SetVarName(2,"prob");
454 fContAllChargesV0Cv3->SetVarName(3,"#Psi");
455 fContAllChargesV0Cv3->SetVarName(4,"PIDmask");
0f25ad32 456 fContAllChargesV0Cv3->SetVarName(5,"DCAbin");
243fbce7 457 if(fV3) fContAllChargesV0Cv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
458 if(fV3) fContAllChargesV0Cv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
459 if(fV3) fContAllChargesV0Cv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
460 if(fV3) fContAllChargesV0Cv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
461 if(fV3) fContAllChargesV0Cv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
462 if(fV3) fContAllChargesV0Cv3->AddSpecies("d",nPtBinsTOF,binsPtTOF);
463 if(fV3) fContAllChargesV0Cv3->AddSpecies("t",nPtBinsTOF,binsPtTOF);
464 if(fV3) fContAllChargesV0Cv3->AddSpecies("he3",nPtBinsTOF,binsPtTOF);
465 if(fV3) fContAllChargesV0Cv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
a640d925 466 if(fV3) fContAllChargesV0Cv3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
467 if(fV3) fContAllChargesV0Cv3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
468 if(fV3) fContAllChargesV0Cv3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
a31d07c8 469 if(fV3) fContAllChargesV0Cv3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
afa8df58 470
471 fList2->Add(fContAllChargesV0Av3);
472 fList2->Add(fContAllChargesV0Cv3);
473
08371bdc 474 if(fIsMC && fV3){
475 fContAllChargesMCAv3 = new AliFlowVZEROResults("v3mcA",5,binsTOFmcPureMC);
a640d925 476 fContAllChargesMCAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
08371bdc 477 fContAllChargesMCAv3->SetVarRange(1,-1.5,1.5); // charge
478 fContAllChargesMCAv3->SetVarRange(2,0.6,1.0001);// prob
479 fContAllChargesMCAv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
480 fContAllChargesMCAv3->SetVarRange(4,-0.5,1.5); // pid mask
481 fContAllChargesMCAv3->SetVarName(0,"centrality");
482 fContAllChargesMCAv3->SetVarName(1,"charge");
483 fContAllChargesMCAv3->SetVarName(2,"prob");
484 fContAllChargesMCAv3->SetVarName(3,"#Psi");
485 fContAllChargesMCAv3->SetVarName(4,"PIDmask");
486 fContAllChargesMCAv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
487 fContAllChargesMCAv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
488 fContAllChargesMCAv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
489 fContAllChargesMCAv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
490 fContAllChargesMCAv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
491 fContAllChargesMCAv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
492 fList3->Add(fContAllChargesMCAv3);
493
494 fContAllChargesMCCv3 = new AliFlowVZEROResults("v3mcC",5,binsTOFmcPureMC);
a640d925 495 fContAllChargesMCCv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
08371bdc 496 fContAllChargesMCCv3->SetVarRange(1,-1.5,1.5); // charge
497 fContAllChargesMCCv3->SetVarRange(2,0.6,1.0001);// prob
498 fContAllChargesMCCv3->SetVarRange(3,-TMath::Pi()/3,TMath::Pi()/3); // Psi
499 fContAllChargesMCCv3->SetVarRange(4,-0.5,1.5); // pid mask
500 fContAllChargesMCCv3->SetVarName(0,"centrality");
501 fContAllChargesMCCv3->SetVarName(1,"charge");
502 fContAllChargesMCCv3->SetVarName(2,"prob");
503 fContAllChargesMCCv3->SetVarName(3,"#Psi");
504 fContAllChargesMCCv3->SetVarName(4,"PIDmask");
505 fContAllChargesMCCv3->AddSpecies("all",nPtBinsTOF,binsPtTOF);
506 fContAllChargesMCCv3->AddSpecies("pi",nPtBinsTOF,binsPtTOF);
507 fContAllChargesMCCv3->AddSpecies("k",nPtBinsTOF,binsPtTOF);
508 fContAllChargesMCCv3->AddSpecies("pr",nPtBinsTOF,binsPtTOF);
509 fContAllChargesMCCv3->AddSpecies("e",nPtBinsTOF,binsPtTOF);
510 fContAllChargesMCCv3->AddSpecies("mu",nPtBinsTOF,binsPtTOF);
511 fList3->Add(fContAllChargesMCCv3);
512 }
513
afa8df58 514 // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
515 // v2
a640d925 516 fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
517 fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
518 fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
afa8df58 519
520 fList->Add(fHResTPCv0A2);
521 fList->Add(fHResTPCv0C2);
522 fList->Add(fHResv0Cv0A2);
523
524 // v3
a640d925 525 fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
526 fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
527 fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
afa8df58 528
529 fList2->Add(fHResTPCv0A3);
530 fList2->Add(fHResTPCv0C3);
531 fList2->Add(fHResv0Cv0A3);
532
08371bdc 533 // MC as in the dataEP resolution (but using MC tracks)
534 if(fIsMC && fV3){
a640d925 535 fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
536 fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
537 fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
08371bdc 538 fList3->Add(fHResMA2);
539 fList3->Add(fHResMC2);
540 fList3->Add(fHResAC2);
541 }
542 if(fIsMC && fV3){
a640d925 543 fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
544 fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
545 fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
08371bdc 546 fList3->Add(fHResMA3);
547 fList3->Add(fHResMC3);
548 fList3->Add(fHResAC3);
549 }
550
551
afa8df58 552 // V0A and V0C event plane distributions
553 //v2
a640d925 554 fPhiRPv0A = new TH2F("fPhiRPv0Av2","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
555 fPhiRPv0C = new TH2F("fPhiRPv0Cv2","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/2,TMath::Pi()/2);
afa8df58 556
557 //v3
a640d925 558 fPhiRPv0Av3 = new TH2F("fPhiRPv0Av3","#phi distribution of EP VZERO-A;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
559 fPhiRPv0Cv3 = new TH2F("fPhiRPv0Cv3","#phi distribution of EP VZERO-C;centrality;#phi (rad)",nCentrBin,0,nCentrBin,nPsiTOF,-TMath::Pi()/3,TMath::Pi()/3);
afa8df58 560
561 // QA container
562 // v2
563 const Int_t nDETsignal = 50;
564 Double_t binDETsignal[nDETsignal+1];
565 for(Int_t i=0;i<nDETsignal+1;i++){
566 binDETsignal[i] = -5 + i*10. / nDETsignal;
567 }
587d006a 568// const Int_t nEta = 5;
569// Double_t binEta[nEta+1];
570// for(Int_t i=0;i<nEta+1;i++){
571// binEta[i] = -1 + i*2. / nEta;
572// }
afa8df58 573
574 const Int_t nDeltaPhi = 5;
afa8df58 575 const Int_t nDeltaPhiV3 = 7;
afa8df58 576
243fbce7 577 Int_t binsQA[5] = {nCentrTOF,7,5,nDeltaPhi,2};
578 Int_t binsQAv3[5] = {nCentrTOF,7,5,nDeltaPhiV3,2};
579
580
581 fQA = new AliFlowVZEROQA("v2AQA",5,binsQA);
a640d925 582 fQA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
243fbce7 583 fQA->SetVarRange(1,0,7); // pt
584 fQA->SetVarRange(2,0.,1.0001);// prob
585 fQA->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
586 fQA->SetVarRange(4,-0.5,1.5); // pid mask
587 fQA->SetVarName(0,"centrality");
588 fQA->SetVarName(1,"p_{t}");
589 fQA->SetVarName(2,"prob");
590 fQA->SetVarName(3,"#Psi");
591 fQA->SetVarName(4,"PIDmask");
592 if(fQAsw && fV2) fQA->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
593 if(fQAsw && fV2) fQA->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
594 if(fQAsw && fV2) fQA->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
595// if(fQAsw && fV2) fQA->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
596// fQA->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
597// fQA->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
598// fQA->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
599
600 fQA2 = new AliFlowVZEROQA("v2CQA",5,binsQA);
a640d925 601 fQA2->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
243fbce7 602 fQA2->SetVarRange(1,0,7); // pt
603 fQA2->SetVarRange(2,0.,1.0001);// prob
604 fQA2->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
605 fQA2->SetVarRange(4,-0.5,1.5); // pid mask
606 fQA2->SetVarName(0,"centrality");
607 fQA2->SetVarName(1,"p_{t}");
608 fQA2->SetVarName(2,"prob");
609 fQA2->SetVarName(3,"#Psi");
610 fQA2->SetVarName(4,"PIDmask");
611 if(fQAsw && fV2) fQA2->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
612 if(fQAsw && fV2) fQA2->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
613 if(fQAsw && fV2) fQA2->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
614// if(fQAsw && fV2) fQA2->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
615// fQA2->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
616// fQA2->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
617// fQA2->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
618
619 fQAv3 = new AliFlowVZEROQA("v3AQA",5,binsQAv3);
a640d925 620 fQAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
243fbce7 621 fQAv3->SetVarRange(1,0,7); // pt
622 fQAv3->SetVarRange(2,0.,1.0001);// prob
623 fQAv3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
624 fQAv3->SetVarRange(4,-0.5,1.5); // pid mask
625 fQAv3->SetVarName(0,"centrality");
626 fQAv3->SetVarName(1,"p_{t}");
627 fQAv3->SetVarName(2,"prob");
628 fQAv3->SetVarName(3,"#Psi");
629 fQAv3->SetVarName(4,"PIDmask");
630 if(fQAsw && fV3) fQAv3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
631// if(fQAsw && fV3) fQAv3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
632// if(fQAsw && fV3) fQAv3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
633// if(fQAsw && fV2) fQAv3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
634// fQAv3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
635// fQAv3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
636// fQAv3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
637
638 fQA2v3 = new AliFlowVZEROQA("v3CQA",5,binsQAv3);
a640d925 639 fQA2v3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
243fbce7 640 fQA2v3->SetVarRange(1,0,7); // pt
641 fQA2v3->SetVarRange(2,0.,1.0001);// prob
642 fQA2v3->SetVarRange(3,-TMath::Pi(),TMath::Pi()); // Psi
643 fQA2v3->SetVarRange(4,-0.5,1.5); // pid mask
644 fQA2v3->SetVarName(0,"centrality");
645 fQA2v3->SetVarName(1,"p_{t}");
646 fQA2v3->SetVarName(2,"prob");
647 fQA2v3->SetVarName(3,"#Psi");
648 fQA2v3->SetVarName(4,"PIDmask");
649 if(fQAsw && fV3) fQA2v3->AddSpecies("pi",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
650// if(fQAsw && fV3) fQA2v3->AddSpecies("k",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
651// if(fQAsw && fV3) fQA2v3->AddSpecies("pr",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
652// if(fQAsw && fV2) fQA2v3->AddSpecies("e",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
653// fQA2v3->AddSpecies("d",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
654// fQA2v3->AddSpecies("t",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
655// fQA2v3->AddSpecies("he3",nDETsignal,binDETsignal,nDETsignal,binDETsignal);
afa8df58 656
657 fList->Add(fPhiRPv0A);
658 fList->Add(fPhiRPv0C);
587d006a 659
660 if(fQAsw && fV2){
661 fList4->Add(fQA);
662 fList4->Add(fQA2);
663 }
afa8df58 664
665 fList2->Add(fPhiRPv0Av3);
666 fList2->Add(fPhiRPv0Cv3);
afa8df58 667
587d006a 668 if(fQAsw && fV3){
669 fList4->Add(fQAv3);
670 fList4->Add(fQA2v3);
671 }
672
a72e0394 673 // fList->Add(fTree); // comment if not needed
afa8df58 674
0f25ad32 675 const Int_t nDCA = 300;
676 Double_t DCAbin[nDCA+1];
677 for(Int_t i=0;i <= nDCA;i++){
678 DCAbin[i] = -3 +i*6.0/nDCA;
679 }
680
681 char nameHistos[100];
682 for(Int_t iC=0;iC < nCentrBin;iC++){
683 snprintf(nameHistos,100,"fHdcaPtPiCent%i",iC);
684 fHdcaPt[iC][0] = new TH2D(nameHistos,"DCA_{xy} for #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
685 snprintf(nameHistos,100,"fHdcaPtKaCent%i",iC);
686 fHdcaPt[iC][1] = new TH2D(nameHistos,"DCA_{xy} for K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
687 snprintf(nameHistos,100,"fHdcaPtPrCent%i",iC);
688 fHdcaPt[iC][2] = new TH2D(nameHistos,"DCA_{xy} for #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
689 snprintf(nameHistos,100,"fHdcaPtElCent%i",iC);
690 fHdcaPt[iC][3] = new TH2D(nameHistos,"DCA_{xy} for e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
691 snprintf(nameHistos,100,"fHdcaPtDeCent%i",iC);
692 fHdcaPt[iC][4] = new TH2D(nameHistos,"DCA_{xy} for #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
693 snprintf(nameHistos,100,"fHdcaPtTrCent%i",iC);
694 fHdcaPt[iC][5] = new TH2D(nameHistos,"DCA_{xy} for #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
695 snprintf(nameHistos,100,"fHdcaPtHeCent%i",iC);
696 fHdcaPt[iC][6] = new TH2D(nameHistos,"DCA_{xy} for #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
697 }
698
699 if(fFillDCA && fQAsw){
700 for(Int_t i=0;i<7;i++)
701 for(Int_t iC=0;iC < nCentrBin;iC++)
702 fList4->Add(fHdcaPt[iC][i]);
703 }
704 if(fIsMC){
705 for(Int_t iC=0;iC < nCentrBin;iC++){
706 snprintf(nameHistos,100,"fHdcaPtPiSecCent%i",iC);
707 fHdcaPtSec[iC][0] = new TH2D(nameHistos,"DCA_{xy} for secondary #pi;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
708 snprintf(nameHistos,100,"fHdcaPtKaSecCent%i",iC);
709 fHdcaPtSec[iC][1] = new TH2D(nameHistos,"DCA_{xy} for secondary K;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
710 snprintf(nameHistos,100,"fHdcaPtPrSecCent%i",iC);
711 fHdcaPtSec[iC][2] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{p};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
712 snprintf(nameHistos,100,"fHdcaPtElSecCent%i",iC);
713 fHdcaPtSec[iC][3] = new TH2D(nameHistos,"DCA_{xy} for secondary e;p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
714 snprintf(nameHistos,100,"fHdcaPtDeSecCent%i",iC);
715 fHdcaPtSec[iC][4] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{d};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
716 snprintf(nameHistos,100,"fHdcaPtTrSecCent%i",iC);
717 fHdcaPtSec[iC][5] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{t};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
718 snprintf(nameHistos,100,"fHdcaPtHeSecCent%i",iC);
719 fHdcaPtSec[iC][6] = new TH2D(nameHistos,"DCA_{xy} for secondary #bar{^{3}He};p_{t} (GeV/c);DCA_{xy} (cm)",nPtBinsTOF,binsPtTOF,nDCA,DCAbin);
720 }
721
722 if(fFillDCA && fQAsw){
723 for(Int_t i=0;i<7;i++)
724 for(Int_t iC=0;iC < nCentrBin;iC++)
725 fList4->Add(fHdcaPtSec[iC][i]);
726 }
727 }
728
729 // Add TProfile Extra QA
730 const Int_t nBinQApid = 2;
731 Int_t binQApid[nBinQApid] = {nCentrTOF,200};
732 const Int_t nbinsigma = 100;
733 Double_t nsigmaQA[nbinsigma+1];
734 for(Int_t i=0;i<nbinsigma+1;i++){
735 nsigmaQA[i] = -10 + 20.0*i/nbinsigma;
736 }
737 fContQApid = new AliFlowVZEROResults("qaPID",nBinQApid,binQApid);
a640d925 738 fContQApid->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
0f25ad32 739 fContQApid->SetVarRange(1,0,20); // charge
740 fContQApid->SetVarName(0,"centrality");
741 fContQApid->SetVarName(1,"p_{t}");
742 fContQApid->AddSpecies("piTPC",nbinsigma,nsigmaQA);
743 fContQApid->AddSpecies("piTOF",nbinsigma,nsigmaQA);
744 fContQApid->AddSpecies("kaTPC",nbinsigma,nsigmaQA);
745 fContQApid->AddSpecies("kaTOF",nbinsigma,nsigmaQA);
746 fContQApid->AddSpecies("prTPC",nbinsigma,nsigmaQA);
747 fContQApid->AddSpecies("prTOF",nbinsigma,nsigmaQA);
748 if(fV2) fList->Add(fContQApid);
749
afa8df58 750 printf("Output creation ok!!\n\n\n\n");
751
a31d07c8 752 fList->Add(fHctauPtEP);
753 fList->Add(fHctauAt1EP);
754
755 fHKsPhi = new TH2D("hKsPhi","K^{0}_{s} #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
756 fList->Add(fHKsPhi);
757 fHKsPhiEP = new TH2D("hKsPhiEP","EP V0C #phi distributuion;v_{z} (cm);#phi (rad)",20,-10,10,20,0,2*TMath::Pi());
758 fList->Add(fHKsPhiEP);
759
760 fHK0sMass = new TH2D("hK0sMass","K^{0}_{s} mass;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
761 fList->Add(fHK0sMass);
762 fHK0sMass2 = new TH2D("hK0sMass2","K^{0}_{s} mass using secondary vertex;p_{T} (GeV/#it{c});mass (GeV/#it{c}^{2})",20,0,5,400,0,1);
763 fList->Add(fHK0sMass2);
764
765 fHK0vsLambda= new TH2D("hK0vsLambda",";K^{0} mass;#Lambda mass",100,0,1,100,0.5,1.5);
766 fList->Add(fHK0vsLambda);
767
afa8df58 768 // Post output data.
769 if(fV2) PostData(1, fList);
770 if(fV3) PostData(2, fList2);
587d006a 771 if(fIsMC) PostData(3, fList3);
772 if(fQAsw) PostData(4, fList4);
243fbce7 773
afa8df58 774}
775
776//______________________________________________________________________________
777void AliAnalysisTaskVnV0::UserExec(Option_t *)
778{
779 // Main loop
780 // Called for each event
781
6fa5bb6e 782 fgIsPsiComputed = kFALSE;
783 fgPsi2v0a=999.;
784 fgPsi2v0c=999.;
785 fgPsi2tpc=999.;
786 fgPsi3v0a=999.;
787 fgPsi3v0c=999.;
788 fgPsi3tpc=999.;
a640d925 789 fgPsi2v0aMC=999.;
790 fgPsi2v0cMC=999.;
791 fgPsi2tpcMC=999.;
792 fgPsi3v0aMC=999.;
793 fgPsi3v0cMC=999.;
794 fgPsi3tpcMC=999.;
6fa5bb6e 795
243fbce7 796 fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
797 if(!fOutputAOD){
afa8df58 798 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
799 this->Dump();
800 return;
801 }
802
243fbce7 803 Int_t run = fOutputAOD->GetRunNumber();
afa8df58 804
805 if(run != fRun){
806 // Load the calibrations run dependent
a31d07c8 807 if(! fIsAfter2011) OpenInfoCalbration(run);
808 fRun=run;
afa8df58 809 }
810
a31d07c8 811 fZvtx = GetVertex(fOutputAOD);
afa8df58 812
587d006a 813
814
815 //Get the MC object
d6a1c304 816 if(fIsMC){
243fbce7 817 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
d6a1c304 818 if (!mcHeader) {
819 AliError("Could not find MC Header in AOD");
820 return;
821 }
587d006a 822 }
823
824 /*
825 AliMCEvent* mcEvent = MCEvent();
826 if (!mcEvent) {
827 Printf("ERROR: Could not retrieve MC event");
828 return;
829 }
830
831 Double_t gReactionPlane = -999., gImpactParameter = -999.;
832 //Get the MC header
833 AliGenHijingEventHeader* headerH = dynamic_cast<AliGenHijingEventHeader*>(mcEvent->GenEventHeader());
834 if (headerH) {
835 //Printf("=====================================================");
836 //Printf("Reaction plane angle: %lf",headerH->ReactionPlaneAngle());
837 //Printf("=====================================================");
838 gReactionPlane = headerH->ReactionPlaneAngle();
839 gImpactParameter = headerH->ImpactParameter();
840 }
841
842*/
843
a31d07c8 844 if (TMath::Abs(fZvtx) < fVtxCut) {
afa8df58 845 //Centrality
846 Float_t v0Centr = -10.;
847 Float_t trkCentr = -10.;
243fbce7 848 AliCentrality *centrality = fOutputAOD->GetCentrality();
afa8df58 849 if (centrality){
a640d925 850// printf("v0centr = %f -- tpccnetr%f\n",centrality->GetCentralityPercentile("V0M"),centrality->GetCentralityPercentile("TRK"));
a31d07c8 851 v0Centr = centrality->GetCentralityPercentile("V0M");
852 trkCentr = centrality->GetCentralityPercentile("TRK");
853 //changed
afa8df58 854 }
855
a31d07c8 856 if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr > 0){ // consistency cut on centrality selection
243fbce7 857 fPID->SetDetResponse(fOutputAOD, v0Centr); // Set the PID object for each event!!!!
858 Analyze(fOutputAOD,v0Centr); // Do analysis!!!!
afa8df58 859
860 fCentrality = v0Centr;
861 if(fV2) fTree->Fill();
862 }
863 }
864
865}
866
867//________________________________________________________________________
868void AliAnalysisTaskVnV0::Analyze(AliAODEvent* aodEvent, Float_t v0Centr)
869{
a31d07c8 870
871 Int_t nusedForK0s=0;
872 AliAODTrack *usedForK0s1[1000];
873 AliAODTrack *usedForK0s2[1000];
874
afa8df58 875 Float_t mass[8] = {5.10998909999999971e-04, 1.05658000000000002e-01, 1.39570000000000000e-01, 4.93676999999999977e-01, 9.38271999999999995e-01,1.87783699999999998,2.81740199999999996,1.40805449999999999};
876
877 // Event plane resolution for v2
878 Float_t evPlRes[18] = {0.350582,0.505393,0.607845,0.632913,0.592230,0.502489,0.381717,0.249539,0.133180, // V0A vs. centrality
879 0.446480,0.612705,0.712222,0.736200,0.697907,0.610114,0.481009,0.327402,0.182277};// V0C vs. centrality
880
881 Int_t iC = -1;
a31d07c8 882 if (v0Centr < 80){ // analysis only for 0-80% centrality classes
883 // if (v0Centr >0 && v0Centr < 80){ // analysis only for 0-80% centrality classes
884 // changed
6fa5bb6e 885 fgIsPsiComputed = kTRUE;
886
afa8df58 887 // centrality bins
888 if(v0Centr < 5) iC = 0;
889 else if(v0Centr < 10) iC = 1;
890 else if(v0Centr < 20) iC = 2;
891 else if(v0Centr < 30) iC = 3;
892 else if(v0Centr < 40) iC = 4;
893 else if(v0Centr < 50) iC = 5;
894 else if(v0Centr < 60) iC = 6;
895 else if(v0Centr < 70) iC = 7;
896 else iC = 8;
a640d925 897
898 Int_t iCcal = iC;
899
a31d07c8 900/*
a640d925 901 if(nCentrBin==16){
a31d07c8 902 iC = v0Centr/5;
a640d925 903 if(iC >= nCentrBin) iC = nCentrBin-1;
904 }
a31d07c8 905
906 // centrality bins
907 // changed
908 if(v0Centr < 10 + 10./9) iC = 0;
909 else if(v0Centr < 10 + 20./9) iC = 1;
910 else if(v0Centr < 10 + 30./9) iC = 2;
911 else if(v0Centr < 10 + 40./9) iC = 3;
912 else if(v0Centr < 10 + 50./9) iC = 4;
913 else if(v0Centr < 10 + 60./9) iC = 5;
914 else if(v0Centr < 10 + 70./9) iC = 6;
915 else if(v0Centr < 10 + 90./9) iC = 7;
916 else if(v0Centr < 10 + 100./9) iC = 8;
917 else if(v0Centr < 10 + 110./9) iC = 9;
918 else if(v0Centr < 10 + 120./9) iC = 10;
919 else if(v0Centr < 10 + 130./9) iC = 11;
920 else if(v0Centr < 10 + 140./9) iC = 12;
921 else if(v0Centr < 10 + 150./9) iC = 13;
922 else if(v0Centr < 10 + 160./9) iC = 14;
923 else if(v0Centr < 10 + 170./9) iC = 15;
924 else iC = 16;
925 if(iC >= nCentrBin) iC= nCentrBin - 1;
926*/
927
afa8df58 928 //reset Q vector info
929 Double_t Qxa2 = 0, Qya2 = 0;
930 Double_t Qxc2 = 0, Qyc2 = 0;
931 Double_t Qxa3 = 0, Qya3 = 0;
932 Double_t Qxc3 = 0, Qyc3 = 0;
587d006a 933
934 Int_t nAODTracks = aodEvent->GetNumberOfTracks();
935
936 AliAODMCHeader *mcHeader = NULL;
937 TClonesArray *mcArray = NULL;
938 Float_t evplaneMC = 0;
939 if(fIsMC){
243fbce7 940 mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
587d006a 941
942 if (mcHeader) {
4edd7307 943 evplaneMC = mcHeader->GetReactionPlaneAngle();
944 if(evplaneMC > TMath::Pi()/2 && evplaneMC <= TMath::Pi()*3/2) evplaneMC-=TMath::Pi();
945 else if(evplaneMC > TMath::Pi()*3/2) evplaneMC-=2*TMath::Pi();
243fbce7 946 mcArray = (TClonesArray*)fOutputAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
08371bdc 947
948 if(mcArray){
949 Float_t QxMCv2[3] = {0,0,0};
950 Float_t QyMCv2[3] = {0,0,0};
951 Float_t QxMCv3[3] = {0,0,0};
952 Float_t QyMCv3[3] = {0,0,0};
953 Float_t EvPlaneMCV2[3] = {0,0,0};
954 Float_t EvPlaneMCV3[3] = {0,0,0};
955 Float_t etaMin[3] = {2.8,-3.6,-0.8}; // A-side, C-side M-barrel
956 Float_t etaMax[3] = {4.88,-1.8,0.8};
957
958 // analysis on MC tracks
959 Int_t nMCtrack = mcArray->GetEntries() ;
960
961 // EP computation with MC tracks
962 for(Int_t iT=0;iT < nMCtrack;iT++){
963 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
964 if(!mctr || !(mctr->IsPrimary()) || !(mctr->Charge()) || mctr->Pt() < 0.2) continue;
965
966 Float_t eta = mctr->Eta();
967
968 for(Int_t iD=0;iD<3;iD++){
969 if(eta > etaMin[iD] && eta < etaMax[iD]){
970 Float_t phi = mctr->Phi();
971 QxMCv2[iD] += TMath::Cos(2*phi);
972 QyMCv2[iD] += TMath::Sin(2*phi);
973 QxMCv3[iD] += TMath::Cos(3*phi);
974 QyMCv3[iD] += TMath::Sin(3*phi);
975 }
976 }
977 }
978
979 if(fV2){
a640d925 980 EvPlaneMCV2[0] = TMath::ATan2(QyMCv2[0],QxMCv2[0])/2.;
981 EvPlaneMCV2[1] = TMath::ATan2(QyMCv2[1],QxMCv2[1])/2.;
982 EvPlaneMCV2[2] = TMath::ATan2(QyMCv2[2],QxMCv2[2])/2.;
08371bdc 983 fHResMA2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[0])));
984 fHResMC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[2]-EvPlaneMCV2[1])));
985 fHResAC2->Fill(Double_t(iC), TMath::Cos(2*(EvPlaneMCV2[0]-EvPlaneMCV2[1])));
a640d925 986 fgPsi2v0aMC = EvPlaneMCV2[0];
987 fgPsi2v0cMC = EvPlaneMCV2[1];
988 fgPsi2tpcMC = EvPlaneMCV2[2];
08371bdc 989 }
990 if(fV3){
a640d925 991 EvPlaneMCV3[0] = TMath::ATan2(QyMCv3[0],QxMCv3[0])/3.;
992 EvPlaneMCV3[1] = TMath::ATan2(QyMCv3[1],QxMCv3[1])/3.;
993 EvPlaneMCV3[2] = TMath::ATan2(QyMCv3[2],QxMCv3[2])/3.;
08371bdc 994 fHResMA3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[0])));
995 fHResMC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[2]-EvPlaneMCV3[1])));
996 fHResAC3->Fill(Double_t(iC), TMath::Cos(3*(EvPlaneMCV3[0]-EvPlaneMCV3[1])));
a640d925 997 fgPsi3v0aMC = EvPlaneMCV3[0];
998 fgPsi3v0cMC = EvPlaneMCV3[1];
999 fgPsi3tpcMC = EvPlaneMCV3[2];
08371bdc 1000 }
1001
1002 // flow A and C side
a31d07c8 1003 Float_t xMCepAv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[0],1};
1004 Float_t xMCepCv2[5] = {iC,0/*charge*/,1,EvPlaneMCV2[1],1};
1005 Float_t xMCepAv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[0],1};
1006 Float_t xMCepCv3[5] = {iC,0/*charge*/,1,EvPlaneMCV3[1],1};
08371bdc 1007
1008 for(Int_t iT=0;iT < nMCtrack;iT++){
1009 AliAODMCParticle *mctr = (AliAODMCParticle*) mcArray->At(iT);
1010 if(!mctr || !(mctr->IsPhysicalPrimary()) || !(mctr->Charge()) || TMath::Abs(mctr->Eta()) > 0.8 || mctr->Pt() < 0.2) continue;
1011 Int_t iS = TMath::Abs(mctr->GetPdgCode());
1012 Int_t charge = mctr->Charge();
1013 Float_t pt = mctr->Pt();
1014 Float_t phi = mctr->Phi();
1015
1016 if(charge > 0){
1017 xMCepAv2[1] = 1;
1018 xMCepCv2[1] = 1;
1019 xMCepAv3[1] = 1;
1020 xMCepCv3[1] = 1;
1021 }
1022 else{
1023 xMCepAv2[1] = -1;
1024 xMCepCv2[1] = -1;
1025 xMCepAv3[1] = -1;
1026 xMCepCv3[1] = -1;
1027 }
1028
a640d925 1029 fContAllChargesMCA->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1030 fContAllChargesMCC->Fill(0,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1031 fContAllChargesMCAv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1032 fContAllChargesMCCv3->Fill(0,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1033
1034 if(iS==11){
a640d925 1035 fContAllChargesMCA->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1036 fContAllChargesMCC->Fill(4,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1037 fContAllChargesMCAv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1038 fContAllChargesMCCv3->Fill(4,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1039 }
1040 else if(iS==13){
a640d925 1041 fContAllChargesMCA->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1042 fContAllChargesMCC->Fill(5,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1043 fContAllChargesMCAv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1044 fContAllChargesMCCv3->Fill(5,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1045 }
1046 else if(iS==211){
a640d925 1047 fContAllChargesMCA->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1048 fContAllChargesMCC->Fill(1,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1049 fContAllChargesMCAv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1050 fContAllChargesMCCv3->Fill(1,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1051 }
1052 else if(iS==321){
a640d925 1053 fContAllChargesMCA->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1054 fContAllChargesMCC->Fill(2,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1055 fContAllChargesMCAv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1056 fContAllChargesMCCv3->Fill(2,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1057 }
1058 else if(iS==2212){
a640d925 1059 fContAllChargesMCA->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[0])),xMCepAv2);
1060 fContAllChargesMCC->Fill(3,pt, TMath::Cos(2*(phi - EvPlaneMCV2[1])),xMCepCv2);
1061 fContAllChargesMCAv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[0])),xMCepAv3);
1062 fContAllChargesMCCv3->Fill(3,pt, TMath::Cos(3*(phi - EvPlaneMCV3[1])),xMCepCv3);
08371bdc 1063 }
1064 }
1065 }
587d006a 1066 }
1067 }
1068
a31d07c8 1069 // TPC EP needed for resolution studies (TPC subevent)
1070 Double_t Qx2 = 0, Qy2 = 0;
1071 Double_t Qx3 = 0, Qy3 = 0;
1072
1073 for(Int_t iT = 0; iT < nAODTracks; iT++) {
1074
1075 AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
1076
1077 if (!aodTrack){
a31d07c8 1078 continue;
1079 }
1080
1081 Bool_t trkFlag = aodTrack->TestFilterBit(1);
1082
1083 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
1084 continue;
1085
1086 Double_t b[2] = {-99., -99.};
1087 Double_t bCov[3] = {-99., -99., -99.};
d8884981 1088
1089
7ae044be 1090 AliAODTrack param(*aodTrack);
1091 if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
a31d07c8 1092 continue;
d8884981 1093 }
a31d07c8 1094
1095 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
1096 continue;
1097
1098 Qx2 += TMath::Cos(2*aodTrack->Phi());
1099 Qy2 += TMath::Sin(2*aodTrack->Phi());
1100 Qx3 += TMath::Cos(3*aodTrack->Phi());
1101 Qy3 += TMath::Sin(3*aodTrack->Phi());
1102
1103 }
1104
1105 evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
1106 evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
1107
1108 fgPsi2tpc = evPlAng2;
1109 fgPsi3tpc = evPlAng3;
1110
1111 SelectK0s();
1112
afa8df58 1113 //V0 info
1114 AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1115
1116 for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1117 Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1118 Float_t multv0 = aodV0->GetMultiplicity(iv0);
a31d07c8 1119
1120 if(! fIsAfter2011){
1121 if(! fIsMC){
1122 if (iv0 < 32){ // V0C
1123 Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1124 Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1125 Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1126 Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1127 } else { // V0A
1128 Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1129 Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1130 Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1131 Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1132 }
1133 }
1134 else{
1135 if (iv0 < 32){ // V0C
1136 Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1137 Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1138 Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1139 Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1140 } else { // V0A
1141 Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1142 Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1143 Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1144 Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1145 }
1146 }
afa8df58 1147 }
1148 }
1149
1150 //grab for each centrality the proper histo with the Qx and Qy to do the recentering
a640d925 1151 Double_t Qxamean2 = fMeanQ[iCcal][1][0];
1152 Double_t Qxarms2 = fWidthQ[iCcal][1][0];
1153 Double_t Qyamean2 = fMeanQ[iCcal][1][1];
1154 Double_t Qyarms2 = fWidthQ[iCcal][1][1];
1155 Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
1156 Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
1157 Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
1158 Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
afa8df58 1159
a640d925 1160 Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
1161 Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
1162 Double_t Qycmean2 = fMeanQ[iCcal][0][1];
1163 Double_t Qycrms2 = fWidthQ[iCcal][0][1];
1164 Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
1165 Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
1166 Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
1167 Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
afa8df58 1168
1169 Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1170 Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1171 Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1172 Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1173 Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
1174 Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
1175 Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
1176 Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
1177
a31d07c8 1178 if(! fIsAfter2011){
1179 if(! fIsMC){
1180 evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
1181 evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
1182 evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
1183 evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
1184 }
1185 else{
1186 evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
1187 evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
1188 evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
1189 evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
1190 }
1191 }
1192 else{
1193 AliEventplane *ep = aodEvent->GetEventplane();
1194 evPlAngV0ACor2 = ep->GetEventplane("V0A", aodEvent, 2);
1195 evPlAngV0CCor2 = ep->GetEventplane("V0C", aodEvent, 2);
1196 evPlAngV0ACor3 = ep->GetEventplane("V0A", aodEvent, 3);
1197 evPlAngV0CCor3 = ep->GetEventplane("V0C", aodEvent, 3);
1198 }
1199
6fa5bb6e 1200
1201 fgPsi2v0a = evPlAngV0ACor2;
6baa5ac2 1202 fgPsi2v0c = evPlAngV0CCor2;
6fa5bb6e 1203 fgPsi3v0a = evPlAngV0ACor3;
6baa5ac2 1204 fgPsi3v0c = evPlAngV0CCor3;
a31d07c8 1205
1206 fHKsPhiEP->Fill(fZvtx,fgPsi2v0c);
afa8df58 1207 //loop track and get pid
1208 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1209 AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
1210
1211 if (!aodTrack){
afa8df58 1212 continue;
1213 }
1214
1215 Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
a31d07c8 1216 if(fFillDCA)
1217 trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
afa8df58 1218
a640d925 1219 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
afa8df58 1220 continue;
1221 }
1222
1223 Double_t b[2] = {-99., -99.};
1224 Double_t bCov[3] = {-99., -99., -99.};
d8884981 1225
1226 AliAODTrack *param = new AliAODTrack(*aodTrack);
1227 if (!param->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
1228 if(param) delete param;
afa8df58 1229 continue;
d8884981 1230 }
1231 if(param) delete param;
afa8df58 1232
0f25ad32 1233 if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
afa8df58 1234 continue;
a31d07c8 1235
1236 if(fFillDCA && (TMath::Abs(b[0]) > 3.0 || TMath::Abs(b[1]) > 3))
1237 continue;
0f25ad32 1238
afa8df58 1239 // re-map the container in an array to do the analysis for V0A and V0C within a loop
1240 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
587d006a 1241 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
243fbce7 1242 AliFlowVZEROQA *QA[2] = {fQA,fQA2};
afa8df58 1243
1244 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
587d006a 1245 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
243fbce7 1246 AliFlowVZEROQA *QAv3[2] = {fQAv3,fQA2v3};
afa8df58 1247
587d006a 1248 // Fill MC results
1249 if(fIsMC && mcArray){
243fbce7 1250 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
587d006a 1251 Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1252
a31d07c8 1253 Float_t xMC[5] = {iC,aodTrack->Charge(),1,evplaneMC,fPID->GetCurrentMask(1)&&tofMismProbMC < 0.5}; // to fill analysis v2 container
587d006a 1254
587d006a 1255 Float_t v2mc = TMath::Cos(2*(aodTrack->Phi() - evplaneMC));
1256
a640d925 1257 fContAllChargesMC->Fill(0,aodTrack->Pt(),v2mc,xMC);
1258
1259 Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->GetPdgCode());
1260 if(iS==11){
1261 fContAllChargesMC->Fill(4,aodTrack->Pt(),v2mc,xMC);
1262 }
1263 else if(iS==13){
1264 fContAllChargesMC->Fill(5,aodTrack->Pt(),v2mc,xMC);
1265 }
1266 else if(iS==211){
1267 fContAllChargesMC->Fill(1,aodTrack->Pt(),v2mc,xMC);
1268 }
1269 else if(iS==321){
1270 fContAllChargesMC->Fill(2,aodTrack->Pt(),v2mc,xMC);
1271 }
1272 else if(iS==2212){
1273 fContAllChargesMC->Fill(3,aodTrack->Pt(),v2mc,xMC);
587d006a 1274 }
1275 }
1276
afa8df58 1277 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1278
0f25ad32 1279 if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
afa8df58 1280
587d006a 1281 Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1282 Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
a31d07c8 1283
243fbce7 1284 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
587d006a 1285 Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
afa8df58 1286 Float_t *probRead = fPID->GetProb();
1287 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1288 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
a31d07c8 1289 Float_t x[6] = {iC,aodTrack->Charge(),1,evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
1290 Float_t x3[6] = {iC,aodTrack->Charge(),1,evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
0f25ad32 1291
1292 // in case fill DCA info
1293 if(fFillDCA){
1294 if(TMath::Abs(b[0]) > 0.1){
1295 x[5] = 1;
1296 x3[5] = 1;
1297 }
1298 if(TMath::Abs(b[0]) > 0.3){
1299 x[5] = 2;
1300 x3[5] = 2;
1301 }
1302 if(fIsMC && mcArray){
1303 if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
1304 x[5] += 3;
1305 x3[5] += 3;
1306 }
1307 }
1308 }
a640d925 1309
afa8df58 1310 // Fill no PID
587d006a 1311 if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
1312 if(fV3) contV0v3[iV0]->Fill(0,aodTrack->Pt(),v3V0,x3);
1313
afa8df58 1314
afa8df58 1315 Double_t dedxExp[8];
1316 Float_t tof = -1;
1317 Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1318 Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
78a615f6 1319
1320 Float_t nsigmaTPC[8];
1321 Float_t nsigmaTOF[8];
1322
afa8df58 1323 if(aodTrack->GetDetPid()){ // check the PID object is available
78a615f6 1324 for(Int_t iS=0;iS < 8;iS++){
afa8df58 1325 dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS);
78a615f6 1326 nsigmaTPC[iS] = (dedx - dedxExp[iS])/(dedxExp[iS]*0.07);
1327 // printf("TPC %i = %f (%f %f)\n",iS, nsigmaTPC[iS],dedx, dedxExp[iS]);
1328 }
587d006a 1329
afa8df58 1330 if(fPID->GetCurrentMask(1)){ // if TOF is present
1331 Float_t ptrack = aodTrack->P();
1332 tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack);
1333 aodTrack->GetIntegratedTimes(inttimes);
1334
1335 for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei
1336 inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]);
1337
78a615f6 1338 for(Int_t iS=0;iS<8;iS++){
1339 expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]);
1340 nsigmaTOF[iS] = (tof - inttimes[iS])/expTOFsigma[iS];
1341 // printf("TOF %i = %f\n",iS, nsigmaTOF[iS]);
1342 }
afa8df58 1343 }
1344 }
1345
1346 Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0];
1347 if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1348 else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1349 if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1350 else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1351
1352 Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0];
1353 if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1354 else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1355 if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1356 else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1357
1358 // variable to fill QA container
a31d07c8 1359 Float_t xQA[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0,x[4]}; // v2
1360 Float_t xQA3[5] = {iC,aodTrack->Pt(), 0.0,deltaPhiV0v3,x[4]}; // v3
afa8df58 1361
0f25ad32 1362 // extra QA TProfiles
1363 if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
a31d07c8 1364 Float_t xQApid[2] = {iC,aodTrack->Pt()};
0f25ad32 1365 fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
1366 fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
1367 fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
1368 fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
1369 fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
1370 fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
1371 }
1372
78a615f6 1373 // QA fill
243fbce7 1374 if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
a640d925 1375 else{
243fbce7 1376 if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
1377 xQA[2] = prob[2];
1378 xQA3[2] = xQA[2];
a31d07c8 1379 if(fQAsw && fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
1380 if(fQAsw && fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
78a615f6 1381 }
243fbce7 1382 if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
1383 xQA[2] = prob[3];
1384 xQA3[2] = xQA[2];
a31d07c8 1385 if(fQAsw && fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
1386// if(fQAsw && fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);
78a615f6 1387 }
243fbce7 1388 if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
1389 xQA[2] = prob[4];
1390 xQA3[2] = xQA[2];
a31d07c8 1391 if(fQAsw && fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
1392// if(fQAsw && fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);
78a615f6 1393 }
243fbce7 1394 if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
1395 xQA[2] = prob[0];
1396 xQA3[2] = xQA[2];
a31d07c8 1397// if(fQAsw && fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
1398// if(fQAsw && fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);
78a615f6 1399 }
243fbce7 1400 if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
1401 xQA[2] = prob[5];
1402 xQA3[2] = xQA[2];
a31d07c8 1403 // if(fQAsw && fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
1404 // if(fQAsw && fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);
78a615f6 1405 }
243fbce7 1406 if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
1407 xQA[2] = prob[6];
1408 xQA3[2] = xQA[2];
a31d07c8 1409 // if(fQAsw && fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
1410 // if(fQAsw && fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);
78a615f6 1411 }
243fbce7 1412 if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
1413 xQA[2] = prob[7];
1414 xQA3[2] = xQA[2];
a31d07c8 1415 // if(fQAsw && fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
1416 // if(fQAsw && fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);
78a615f6 1417 }
78a615f6 1418 }
1419
afa8df58 1420 //pid selection
1421 if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!)
1422 else if(prob[2] > 0.6){ // pi
587d006a 1423 x[2] = prob[2];
587d006a 1424 x3[2] = x[2];
78a615f6 1425 if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
587d006a 1426 if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
a72e0394 1427 if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
0f25ad32 1428 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
1429 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1430 }
1431 }
1432 else if(prob[3] > 0.6){ // K
587d006a 1433 x[2] = prob[3];
587d006a 1434 x3[2] = x[2];
78a615f6 1435 if(TMath::Abs(nsigmaTPC[3]) < 5){
587d006a 1436 if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
a72e0394 1437 if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
0f25ad32 1438 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
1439 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1440 }
1441 }
1442 else if(prob[4] > 0.6){ // p
587d006a 1443 x[2] = prob[4];
587d006a 1444 x3[2] = x[2];
78a615f6 1445 if(TMath::Abs(nsigmaTPC[4]) < 5){
587d006a 1446 if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
a72e0394 1447 if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
0f25ad32 1448 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
1449 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1450 }
1451 }
1452 else if(prob[0] > 0.6){ // e
587d006a 1453 x[2] = prob[0];
587d006a 1454 x3[2] = x[2];
78a615f6 1455 if(TMath::Abs(nsigmaTPC[0]) < 5){
587d006a 1456 if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
a72e0394 1457 if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
0f25ad32 1458 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
1459 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1460 }
1461 }
587d006a 1462 else if(prob[1] > 0.6){ // mu
587d006a 1463 x[2] = prob[1];
587d006a 1464 x3[2] = x[2];
78a615f6 1465 if(TMath::Abs(nsigmaTPC[1]) < 5){
587d006a 1466 if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
a72e0394 1467 if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
587d006a 1468 }
1469 }
afa8df58 1470 else if(prob[5] > 0.6){ // d
587d006a 1471 x[2] = prob[5];
587d006a 1472 x3[2] = x[2];
78a615f6 1473 if(TMath::Abs(nsigmaTPC[5]) < 5){
587d006a 1474 if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
a72e0394 1475 if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
0f25ad32 1476 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
1477 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1478 }
1479 }
1480 else if(prob[6] > 0.6){ // t
587d006a 1481 x[2] = prob[6];
587d006a 1482 x3[2] = x[2];
78a615f6 1483 if(TMath::Abs(nsigmaTPC[6]) < 5){
587d006a 1484 if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
a72e0394 1485 if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
0f25ad32 1486 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
1487 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1488 }
1489 }
1490 else if(prob[7] > 0.6){ // He3
587d006a 1491 x[2] = prob[7];
587d006a 1492 x3[2] = x[2];
78a615f6 1493 if(TMath::Abs(nsigmaTPC[7]) < 5){
587d006a 1494 if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
a72e0394 1495 if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
0f25ad32 1496 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
1497 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1498 }
afa8df58 1499 }
1500
587d006a 1501 if(x[4] > 0.5){ // if TOF was present redo TPC stand alone PID to check the PID in the same acceptance (PID mask = 2)
afa8df58 1502 fPID->ResetDetOR(1); // exclude TOF from PID
1503 tofMismProb = 0;
1504
243fbce7 1505 fPID->ComputeProb(aodTrack,fOutputAOD);
587d006a 1506 dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
afa8df58 1507 probRead = fPID->GetProb();
1508
1509 fPID->SetDetOR(1); // include TOF for PID
1510 }
1511 Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities
1512
1513 //pid selection TPC S.A. with TOF matching
587d006a 1514 x[4]*=2; // set the mask to 2 id TOF is present
1515 if(x[4]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance
afa8df58 1516 else if(probTPC[2] > 0.6){ // pi
587d006a 1517 x[2] = probTPC[2];
1518 x3[2] = x[2];
78a615f6 1519 if(TMath::Abs(nsigmaTPC[2]) < 5){
587d006a 1520 if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
a72e0394 1521 if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
afa8df58 1522 }
1523 }
1524 else if(probTPC[3] > 0.6){ // K
587d006a 1525 x[2] = probTPC[3];
1526 x3[2] = x[2];
78a615f6 1527 if(TMath::Abs(nsigmaTPC[3]) < 5){
587d006a 1528 if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
a72e0394 1529 if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
afa8df58 1530 }
1531 }
1532 else if(probTPC[4] > 0.6){ // p
587d006a 1533 x[2] = probTPC[4];
1534 x3[2] = x[2];
78a615f6 1535 if(TMath::Abs(nsigmaTPC[4]) < 5){
587d006a 1536 if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
a72e0394 1537 if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
afa8df58 1538 }
1539 }
1540 else if(probTPC[0] > 0.6){ // e
587d006a 1541 x[2] = probTPC[0];
1542 x3[2] = x[2];
78a615f6 1543 if(TMath::Abs(nsigmaTPC[0]) < 5){
587d006a 1544 if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
a72e0394 1545 if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
587d006a 1546 }
1547 }
1548 else if(probTPC[1] > 0.6){ // mu
1549 x[2] = probTPC[1];
1550 x3[2] = x[2];
78a615f6 1551 if(TMath::Abs(nsigmaTPC[1]) < 5){
587d006a 1552 if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
a72e0394 1553 if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
afa8df58 1554 }
1555 }
1556 else if(probTPC[5] > 0.6){ // d
587d006a 1557 x[2] = probTPC[5];
1558 x3[2] = x[2];
78a615f6 1559 if(TMath::Abs(nsigmaTPC[5]) < 5){
587d006a 1560 if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
a72e0394 1561 if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
afa8df58 1562 }
1563 }
1564 else if(probTPC[6] > 0.6){ // t
587d006a 1565 x[2] = probTPC[6];
1566 x3[2] = x[2];
78a615f6 1567 if(TMath::Abs(nsigmaTPC[6]) < 5){
587d006a 1568 if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
a72e0394 1569 if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
afa8df58 1570 }
1571 }
1572 else if(probTPC[7] > 0.6){ // He3
587d006a 1573 x[2] = probTPC[7];
1574 x3[2] = x[2];
78a615f6 1575 if(TMath::Abs(nsigmaTPC[7]) < 5){
587d006a 1576 if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
a72e0394 1577 if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
afa8df58 1578 }
afa8df58 1579 }
1580 } // end side loop
1581 } // end track loop
1582
a31d07c8 1583 // my V0 loop
1584 for(Int_t imy=0;imy<fNK0s;imy++){
1585 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1586 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1587
1588 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1589 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1590
1591 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1592 Float_t x[6] = {iC,-1/*my K0s are negative for convention*/,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1593 Float_t x3[6] = {iC,-1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
1594
1595 Float_t v2V0 = TMath::Cos(2*(fPhiK0s[imy] - evPlAngV0[iV0]));
1596 Float_t v3V0 = TMath::Cos(3*(fPhiK0s[imy] - evPlAngV0v3[iV0]));
1597 if(fV2) contV0[iV0]->Fill(9,fPtK0s[imy],v2V0,x);
1598 if(fV3) contV0v3[iV0]->Fill(9,fPtK0s[imy],v3V0,x3);
1599 }
1600 }
1601
a640d925 1602 // V0 loop
1603 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1604 AliAODv0 *myV0;
a31d07c8 1605 Double_t dQT, dALPHA, dPT, dMASS=0.0;
a640d925 1606 for (Int_t i=0; i!=nV0s; ++i) {
1607 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1608 if(!myV0) continue;
1609 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > fEtaCut) continue; // skipping low momentum
1610 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0);
1611 if(pass) {
1612 dMASS = myV0->MassK0Short();
1613 pass = 3;
a31d07c8 1614 fHK0sMass2->Fill(myV0->Pt(),dMASS);
a640d925 1615 }
a31d07c8 1616 if(TMath::Abs(dMASS-0.497)/0.005 > 3){
a640d925 1617 pass = PassesAODCuts(myV0,fOutputAOD,1);
1618 if(pass) dMASS = myV0->MassLambda();
1619 if(pass==2) dMASS = myV0->MassAntiLambda();
1620 }
a31d07c8 1621
a640d925 1622 if(pass){// 1 lambda, 2 antilambda, 3=K0s
a31d07c8 1623 dPT=myV0->Pt();
1624 dQT=myV0->PtArmV0();
1625 dALPHA=myV0->AlphaV0();
a640d925 1626
1627 Int_t iPos, iNeg;
1628 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0);
1629 if(iT->Charge()>0) {
1630 iPos = 0; iNeg = 1;
1631 } else {
1632 iPos = 1; iNeg = 0;
1633 }
a31d07c8 1634
1635 // check if one of the daugthers was already used
1636 if(pass == 3 && TMath::Abs(dMASS-0.497)/0.005 < 1){
1637 fHKsPhi->Fill(fZvtx, myV0->Phi());
1638 }
1639
1640 if(pass == 1000){ // disable
1641 Bool_t used = kFALSE;
1642 for(Int_t ii=0;ii<nusedForK0s;ii++){
1643 if(myV0->GetDaughter(iNeg) == usedForK0s1[ii] || myV0->GetDaughter(iPos) == usedForK0s2[ii]){
1644 used = kTRUE;
1645 }
1646 }
1647 if((!used) && nusedForK0s < 1000){
1648 nusedForK0s++;
1649 usedForK0s1[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iNeg);
1650 usedForK0s2[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iPos);
1651 printf("accepted\n");
1652 }
1653 else{
1654 dMASS = 0;
1655 printf("rejected\n");
1656 }
1657 }
1658
a640d925 1659 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1660 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1661
1662 // re-map the container in an array to do the analysis for V0A and V0C within a loop
1663 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1664 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1665
1666 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1667 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1668
1669 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1670
a31d07c8 1671 if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
1672
a640d925 1673 Float_t v2V0 = TMath::Cos(2*(myV0->Phi() - evPlAngV0[iV0]));
1674 Float_t v3V0 = TMath::Cos(3*(myV0->Phi() - evPlAngV0v3[iV0]));
1675
a31d07c8 1676 Float_t deltaphi = myV0->Phi()- evPlAngV0[iV0];
1677 if(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1678 if(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1679
1680 Float_t x[6] = {iC,1,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1681 Float_t x3[6] = {iC,1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
a640d925 1682
1683 Float_t decaylength = myV0->DecayLengthXY(fOutputAOD->GetPrimaryVertex());
1684 // printf("decay length = %f\n",decaylength);
1685
1686 if(pass==2){ // anti-lambda charge = -1
1687 x[1] = -1;
1688 x3[1] = -1;
1689 }
1690
1691 if(decaylength < fMinDistV0) pass = 0;
1692 if(decaylength > fMaxDistV0) pass = 0;
1693
1694 Float_t nsigma = 0;
1695 if(pass < 3)
1696 nsigma = TMath::Abs(dMASS-1.116)/0.0016;
1697 else if(pass == 3)
1698 nsigma = TMath::Abs(dMASS-0.497)/0.005;
1699
1700 if(nsigma < 1)
1701 x[2] = 0.95;
1702 else if(nsigma < 2)
1703 x[2] = 0.85;
1704 else if(nsigma < 3)
1705 x[2] = 0.75;
1706 else if(nsigma < 4)
1707 x[2] = 0.65;
1708 else
1709 x[2] = 0.5;
1710
1711 x3[2] = x[2];
1712
1713 // Fill Container for lambda and Ks
a31d07c8 1714 if(fV2 && pass == 3 && x[2] > 0.6){
1715 contV0[iV0]->Fill(9,myV0->Pt(),v2V0,x);
1716 fHctauPtEP->Fill(myV0->Pt(),deltaphi,decaylength);//ciao
1717 if(myV0->Pt() < 1.1 && myV0->Pt() > 0.9) fHctauAt1EP->Fill(decaylength,deltaphi);
1718 }
a640d925 1719 if(fV3 && pass == 3 && x[2] > 0.6) contV0v3[iV0]->Fill(9,myV0->Pt(),v3V0,x3);
1720 if(fV2 && pass < 3 && x[2] > 0.6) contV0[iV0]->Fill(10,myV0->Pt(),v2V0,x);
1721 if(fV3 && pass < 3 && x[2] > 0.6) contV0v3[iV0]->Fill(10,myV0->Pt(),v3V0,x3);
1722
1723 if(pass < 3){ // lambda
1724 AliAODTrack* aodTrack = iT;
1725 if(pass==2) aodTrack=jT;
1726
1727 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1728 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1729
1730 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1731 Float_t *probRead = fPID->GetProb();
1732 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1733 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
a31d07c8 1734
1735 if(prob[4] < 0.61) prob[4] = 0.61;
a640d925 1736
a31d07c8 1737 Float_t xdec[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v2 container
1738 Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[4],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to fill analysis v3 container
a640d925 1739
1740 // Fill Container for (anti)proton from lambda
a31d07c8 1741 if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
a640d925 1742 if(fV2) contV0[iV0]->Fill(11,aodTrack->Pt(),v2V0,xdec);
1743 if(fV3) contV0v3[iV0]->Fill(11,aodTrack->Pt(),v3V0,xdec3);
1744 }
1745 }
a31d07c8 1746 else if(pass == 3){
1747 AliAODTrack* aodTrack = iT;
1748
1749 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1750 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1751
1752 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1753 Float_t *probRead = fPID->GetProb();
1754 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1755 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1756
1757 if(prob[2] < 0.61) prob[2] = 0.61;
1758
1759 Float_t xdec[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to
1760 Float_t xdec3[6] = {iC,aodTrack->Charge(),prob[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb < 0.5,0}; // to
1761
1762 if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1763 if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdec);
1764 if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdec3);
1765 }
1766
1767 aodTrack = jT;
1768 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1769 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1770
1771 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1772 Float_t *probRead2 = fPID->GetProb();
1773 Float_t prob2[8] = {probRead2[0],probRead2[1],probRead2[2],probRead2[3],probRead2[4],probRead2[5],probRead2[6],probRead2[7]};
1774 Float_t tofMismProb2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1775
1776 if(prob2[2] < 0.61) prob2[2] = 0.61;
1777
1778 Float_t xdecB[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
1779 Float_t xdecB3[6] = {iC,aodTrack->Charge(),prob2[2],evPlAngV0v3[iV0],fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5,0}; // to
1780
1781 if(nsigma < 2 && xdecB[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1782 if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdecB);
1783 if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdecB3);
1784 }
1785 }
a640d925 1786 }
1787
1788 }
1789 } // end loop on V0
1790
1791
afa8df58 1792 // Fill EP distribution histograms
1793 if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
1794 if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
1795
1796 if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
1797 if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
1798
afa8df58 1799 // Fill histograms needed for resolution evaluation
1800 if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
1801 if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
1802 if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
1803
1804 if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
1805 if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
1806 if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
1807 }
1808
a31d07c8 1809
1810
1811 // clean track array
1812 for(Int_t i=0;i < nusedForK0s;i++){
1813 usedForK0s1[i] = NULL;
1814 usedForK0s2[i] = NULL;
1815 }
afa8df58 1816}
1817
1818//_____________________________________________________________________________
1819Float_t AliAnalysisTaskVnV0::GetVertex(AliAODEvent* aod) const
1820{
1821
1822 Float_t zvtx = -999;
1823
1824 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
1825 if (!vtxAOD)
1826 return zvtx;
1827 if(vtxAOD->GetNContributors()>0)
1828 zvtx = vtxAOD->GetZ();
1829
1830 return zvtx;
1831}
1832//_____________________________________________________________________________
1833void AliAnalysisTaskVnV0::Terminate(Option_t *)
1834{
1835 // Terminate loop
1836 Printf("Terminate()");
1837}
1838//_____________________________________________________________________________
1839void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
1840 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1841 TFile *foadb = TFile::Open(oadbfilename.Data());
1842
1843 if(!foadb){
1844 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1845 return;
1846 }
1847
1848 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1849 if(!cont){
1850 printf("OADB object hMultV0BefCorr is not available in the file\n");
1851 return;
1852 }
1853
1854 if(!(cont->GetObject(run))){
1855 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1856 run = 137366;
1857 }
1858 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1859
1860 TF1 *fpol0 = new TF1("fpol0","pol0");
1861 fMultV0->Fit(fpol0,"","",0,31);
1862 fV0Cpol = fpol0->GetParameter(0);
1863 fMultV0->Fit(fpol0,"","",32,64);
1864 fV0Apol = fpol0->GetParameter(0);
1865
1866 for(Int_t iside=0;iside<2;iside++){
1867 for(Int_t icoord=0;icoord<2;icoord++){
a640d925 1868 for(Int_t i=0;i < 9;i++){
afa8df58 1869 char namecont[100];
1870 if(iside==0 && icoord==0)
d6a1c304 1871 snprintf(namecont,100,"hQxc2_%i",i);
afa8df58 1872 else if(iside==1 && icoord==0)
d6a1c304 1873 snprintf(namecont,100,"hQxa2_%i",i);
afa8df58 1874 else if(iside==0 && icoord==1)
d6a1c304 1875 snprintf(namecont,100,"hQyc2_%i",i);
afa8df58 1876 else if(iside==1 && icoord==1)
d6a1c304 1877 snprintf(namecont,100,"hQya2_%i",i);
afa8df58 1878
1879 cont = (AliOADBContainer*) foadb->Get(namecont);
1880 if(!cont){
1881 printf("OADB object %s is not available in the file\n",namecont);
1882 return;
1883 }
1884
1885 if(!(cont->GetObject(run))){
1886 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1887 run = 137366;
1888 }
1889 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1890 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1891
1892 //for v3
1893 if(iside==0 && icoord==0)
d6a1c304 1894 snprintf(namecont,100,"hQxc3_%i",i);
afa8df58 1895 else if(iside==1 && icoord==0)
d6a1c304 1896 snprintf(namecont,100,"hQxa3_%i",i);
afa8df58 1897 else if(iside==0 && icoord==1)
d6a1c304 1898 snprintf(namecont,100,"hQyc3_%i",i);
afa8df58 1899 else if(iside==1 && icoord==1)
d6a1c304 1900 snprintf(namecont,100,"hQya3_%i",i);
afa8df58 1901
1902 cont = (AliOADBContainer*) foadb->Get(namecont);
1903 if(!cont){
1904 printf("OADB object %s is not available in the file\n",namecont);
1905 return;
1906 }
1907
1908 if(!(cont->GetObject(run))){
1909 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1910 run = 137366;
1911 }
1912 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1913 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1914
1915 }
1916 }
1917 }
1918}
a640d925 1919//=======================================================================
1920Int_t AliAnalysisTaskVnV0::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1921{
a31d07c8 1922 Int_t set = 0;
a640d925 1923 Float_t fV0Cuts[9];
1924 // defines cuts to be used
1925 // fV0Cuts[9] dl dca ctp d0 d0d0 qt minEta maxEta PID
1926 switch(set) {
1927 case(0): // No cuts
1928 fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1929 fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1930 fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 0;
1931 break;
1932 case(1): // Tight cuts
1933 fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1934 fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1935 fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 0;
1936 break;
1937 case(2): // Tight cuts + PID
1938 fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1939 fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1940 fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 1;
1941 break;
1942 case(3): // No cuts + PID
1943 fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1944 fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1945 fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 1;
1946 break;
1947 }
1948
1949 // daughter cuts
1950 if(! fCutsDaughter){
1951 fCutsDaughter = new AliESDtrackCuts(Form("daughter_cuts_%s","ESD") );
1952 fCutsDaughter->SetPtRange(0.2,10.0);
1953 fCutsDaughter->SetEtaRange(-0.8, 0.8 );
1954 fCutsDaughter->SetMinNClustersTPC(80);
1955 fCutsDaughter->SetMaxChi2PerClusterTPC(4.0);
1956 fCutsDaughter->SetRequireTPCRefit(kTRUE);
1957 fCutsDaughter->SetAcceptKinkDaughters(kFALSE);
1958 }
1959
1960 if (myV0->GetOnFlyStatus() ) return 0;
1961 //the following is needed in order to evualuate track-quality
1962 AliAODTrack *iT, *jT;
1963 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1964 Double_t pos[3],cov[6];
1965 vV0s->GetXYZ(pos);
1966 vV0s->GetCovarianceMatrix(cov);
1967 const AliESDVertex vESD(pos,cov,100.,100);
1968 // TESTING CHARGE
1969 int iPos, iNeg;
1970 iT=(AliAODTrack*) myV0->GetDaughter(0);
1971 if(iT->Charge()>0) {
1972 iPos = 0; iNeg = 1;
1973 } else {
1974 iPos = 1; iNeg = 0;
1975 }
1976 // END OF TEST
1977
1978 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1979 AliESDtrack ieT( iT );
1980 ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1981 ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1982 ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1983 ieT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1984 if (!fCutsDaughter->IsSelected( &ieT ) ) return 0;
1985
1986 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1987 AliESDtrack jeT( jT );
1988 jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1989 jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1990 jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1991 jeT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1992 if (!fCutsDaughter->IsSelected( &jeT ) ) return 0;
1993
1994 Double_t pvertex[3];
1995 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1996 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1997 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1998 Double_t dDL=myV0->DecayLengthV0( pvertex );
1999 Double_t dDCA=myV0->DcaV0Daughters();
2000 Double_t dCTP=myV0->CosPointingAngle( pvertex );
2001 Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2002 Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2003 Double_t dD0D0=dD0P*dD0M;
2004 Double_t dQT=myV0->PtArmV0();
2005 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
2006 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
2007// Double_t dPT=myV0->Pt();
2008 Double_t dETA=myV0->Eta();
2009 Int_t passes = 1;
2010 if(dDL <fV0Cuts[0]) passes = 0;
2011 if(dDCA >fV0Cuts[1]) passes = 0;
2012 if(dCTP <fV0Cuts[2]) passes = 0;
2013 if(TMath::Abs(dD0P) <fV0Cuts[3]) passes = 0;
2014 if(TMath::Abs(dD0M) <fV0Cuts[3]) passes = 0;
2015 if(dD0D0>fV0Cuts[4]) passes = 0;
2016 if(dETA <fV0Cuts[6]) passes = 0;
2017 if(dETA >fV0Cuts[7]) passes = 0;
2018 if(specie==0) if(dQT<fV0Cuts[5]) passes = 0;
2019 if(specie==1&&passes==1&&dALPHA<0) passes = 2; // antilambda
a31d07c8 2020
2021
2022// if(jT->Pt() < 0.5*myV0->Pt() || iT->Pt() < 0.5*myV0->Pt()) passes = 0;
2023
2024
2025 // additional cut
2026// if(!(iT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2027// if(!(jT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2028
2029// if(!(iT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2030// if(!(jT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2031
2032// if(!(iT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2033// if(!(jT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2034
2035 Bool_t trkFlag = iT->TestFilterBit(1); // TPC only tracks (4,global track)
2036 Bool_t trkFlag2 = jT->TestFilterBit(1); // TPC only tracks (4,global track)
2037
2038 if(!trkFlag) passes = 0;
2039 if(!trkFlag2) passes = 0;
2040
a640d925 2041 if(passes&&fV0Cuts[8]) {
2042
2043 Double_t dedxExp[8];
2044 fPID->ComputeProb(iT,tAOD); // compute Bayesian probabilities
2045 Float_t nsigmaTPC[8];
a31d07c8 2046
2047 Int_t tofMatch=0;
2048 Int_t tofMatch2=0;
2049
a640d925 2050 if(iT->GetDetPid()){ // check the PID object is available
2051 for(Int_t iS=0;iS < 8;iS++){
2052 dedxExp[iS] = fPID->GetExpDeDx(iT,iS);
2053 nsigmaTPC[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2054 }
2055 }
2056 else{
2057 for(Int_t iS=0;iS < 8;iS++)
2058 nsigmaTPC[iS] = 10;
2059 }
2060
a31d07c8 2061
2062 if(fPID->GetCurrentMask(1)) // if TOF is present
2063 tofMatch = 1;
2064
2065// Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2066
2067 Float_t *probRead = fPID->GetProb();
2068 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2069
a640d925 2070 fPID->ComputeProb(jT,tAOD); // compute Bayesian probabilities
2071 Float_t nsigmaTPC2[8];
2072 if(jT->GetDetPid()){ // check the PID object is available
2073 for(Int_t iS=0;iS < 8;iS++){
2074 dedxExp[iS] = fPID->GetExpDeDx(jT,iS);
2075 nsigmaTPC2[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2076 }
2077 }
2078 else{
2079 for(Int_t iS=0;iS < 8;iS++)
2080 nsigmaTPC2[iS] = 10;
2081 }
2082
a31d07c8 2083 if(fPID->GetCurrentMask(1)) // if TOF is present
2084 tofMatch2 = 1;
2085
2086// Float_t tofMismProbMC2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2087
2088 probRead = fPID->GetProb();
2089 Float_t prob2[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2090
a640d925 2091 if(jT->GetTPCNcls() < fNcluster) passes = 0;
2092 else if(iT->GetTPCNcls() < fNcluster) passes = 0;
2093
a31d07c8 2094// if(! (tofMatch && tofMatch2)) passes = 0;
2095
2096 /*
2097 Float_t dMASS = myV0->MassK0Short();
2098 Float_t nsigmaMass = TMath::Abs(dMASS-0.497)/0.005;
2099 if(specie == 0 && TMath::Abs(nsigmaMass) < 1 && myV0->Pt() > 1) printf("candidate i=(pt=%f-phi=%f-tof=%i) j=(pt=%f-phi=%f-tof=%i) \n",iT->Pt(),iT->Phi(),tofMatch,jT->Pt(),jT->Phi(),tofMatch2);
2100 */
2101
a640d925 2102 switch(specie) {
2103 case 0: // K0 PID
a31d07c8 2104 if(0){
2105 if( ((jT->GetTPCmomentum()<15) &&
2106 (TMath::Abs(nsigmaTPC2[2])>3.)) || prob2[2] < 0.9)
2107 passes = 0;
2108 if( ((iT->GetTPCmomentum()<15) &&
2109 (TMath::Abs(nsigmaTPC[2])>3.))|| prob[2] < 0.9 )
2110 passes = 0;
2111 }
a640d925 2112 break;
2113 case 1: // Lambda PID i==pos j ==neg
2114 if(passes==1) {
2115 if( (iT->GetTPCmomentum()<15) &&
2116 (TMath::Abs(nsigmaTPC[4])>3.) )
2117 passes = 0;
2118 if( (jT->GetTPCmomentum()<15) &&
2119 (TMath::Abs(nsigmaTPC2[2])>3.) )
2120 passes = 0;
2121 }
2122 if(passes==2) {
2123 if( (iT->GetTPCmomentum()<15) &&
2124 (TMath::Abs(nsigmaTPC[2])>3.) )
2125 passes = 0;
2126 if( (jT->GetTPCmomentum()<15) &&
2127 (TMath::Abs(nsigmaTPC2[4])>3.) )
2128 passes = 0;
2129 }
2130 break;
2131 }
2132 }
2133 return passes;
2134}
a31d07c8 2135//=======================================================================
2136void AliAnalysisTaskVnV0::SelectK0s(){
2137 fNK0s=0;
2138 fNpiPos=0;
2139 fNpiNeg=0;
2140
2141 if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAng2,1.0); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
2142
2143 // fill pion stacks
2144 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
2145 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
2146 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
2147
2148 if (!aodTrack){
a31d07c8 2149 continue;
2150 }
2151
2152 Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
2153// trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
2154
2155 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
2156 continue;
2157 }
2158
2159 Double_t b[2] = {-99., -99.};
2160 Double_t bCov[3] = {-99., -99., -99.};
d8884981 2161
2162 AliAODTrack *param = new AliAODTrack(*aodTrack);
2163 if (!param->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
2164 if(param) delete param;
a31d07c8 2165 continue;
d8884981 2166 }
2167 if(param) delete param;
a31d07c8 2168
2169 if(TMath::Abs(b[0]) < 0.5/aodTrack->Pt()) continue;
2170
2171 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
2172 Float_t *probRead = fPID->GetProb();
2173 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2174 // Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2175
2176 Int_t charge = aodTrack->Charge();
2177 if(prob[2] > 0.9){
2178 if(charge > 0){
2179 fIPiPos[fNpiPos] = iT;
2180 fNpiPos++;
2181 }
2182 else{
2183 fIPiNeg[fNpiNeg] = iT;
2184 fNpiNeg++;
2185 }
2186 }
2187 }
2188
2189 for(Int_t i=0;i < fNpiPos;i++){
2190 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
2191 AliESDtrack pipE(pip);
2192
2193 for(Int_t j=0;j < fNpiNeg;j++){
2194 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
2195 AliESDtrack pinE(pin);
2196
2197 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
2198
2199 Double_t pPos[3];
2200 Double_t pNeg[3];
2201 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
2202 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
2203
2204 Float_t length = (xp+xn)*0.5;
2205
2206 Float_t pxs = pPos[0] + pNeg[0];
2207 Float_t pys = pPos[1] + pNeg[1];
2208 Float_t pzs = pPos[2] + pNeg[2];
2209 Float_t es = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
2210
2211 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
2212 Float_t phi = TMath::ATan2(pys,pxs);
2213 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
2214
2215 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
2216
2217 if(mindist < 0.2&& length > 1 && length < 25){
2218 fHK0sMass->Fill(pt,mass);
2219
2220 Float_t esL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.938*0.938) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.13957*0.13957);
2221 Float_t esAL = TMath::Sqrt(pPos[0]*pPos[0] + pPos[1]*pPos[1] + pPos[2]*pPos[2] + 0.13957*0.13957) + TMath::Sqrt(pNeg[0]*pNeg[0] + pNeg[1]*pNeg[1] + pNeg[2]*pNeg[2] + 0.938*0.938);
2222
2223 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
2224 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
2225
2226 fHK0vsLambda->Fill(mass,TMath::Min(massaL,massaAL));
2227
2228 if(TMath::Abs(mass-0.497)/0.005 < 1 && massaL > 1.15 && massaAL > 1.15){
2229 fPhiK0s[fNK0s] = phi;
2230 fPtK0s[fNK0s] = pt;
2231 fNK0s++;
2232 }
2233 }
2234 }
2235 }
2236}