]> 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
coverity fix
[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
2942f542 1003 Float_t xMCepAv2[5] = {static_cast<Float_t>(iC),0/*charge*/,1,EvPlaneMCV2[0],1};
1004 Float_t xMCepCv2[5] = {static_cast<Float_t>(iC),0/*charge*/,1,EvPlaneMCV2[1],1};
1005 Float_t xMCepAv3[5] = {static_cast<Float_t>(iC),0/*charge*/,1,EvPlaneMCV3[0],1};
1006 Float_t xMCepCv3[5] = {static_cast<Float_t>(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
2019aa05 1226 AliAODTrack param(*aodTrack);
1227 if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
1228 continue;
d8884981 1229 }
afa8df58 1230
0f25ad32 1231 if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
afa8df58 1232 continue;
a31d07c8 1233
1234 if(fFillDCA && (TMath::Abs(b[0]) > 3.0 || TMath::Abs(b[1]) > 3))
1235 continue;
0f25ad32 1236
afa8df58 1237 // re-map the container in an array to do the analysis for V0A and V0C within a loop
1238 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
587d006a 1239 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
243fbce7 1240 AliFlowVZEROQA *QA[2] = {fQA,fQA2};
afa8df58 1241
1242 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
587d006a 1243 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
243fbce7 1244 AliFlowVZEROQA *QAv3[2] = {fQAv3,fQA2v3};
afa8df58 1245
587d006a 1246 // Fill MC results
1247 if(fIsMC && mcArray){
243fbce7 1248 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
587d006a 1249 Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1250
2942f542 1251 Float_t xMC[5] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),1,static_cast<Float_t>(evplaneMC),static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProbMC < 0.5)}; // to fill analysis v2 container
587d006a 1252
587d006a 1253 Float_t v2mc = TMath::Cos(2*(aodTrack->Phi() - evplaneMC));
1254
a640d925 1255 fContAllChargesMC->Fill(0,aodTrack->Pt(),v2mc,xMC);
1256
1257 Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->GetPdgCode());
1258 if(iS==11){
1259 fContAllChargesMC->Fill(4,aodTrack->Pt(),v2mc,xMC);
1260 }
1261 else if(iS==13){
1262 fContAllChargesMC->Fill(5,aodTrack->Pt(),v2mc,xMC);
1263 }
1264 else if(iS==211){
1265 fContAllChargesMC->Fill(1,aodTrack->Pt(),v2mc,xMC);
1266 }
1267 else if(iS==321){
1268 fContAllChargesMC->Fill(2,aodTrack->Pt(),v2mc,xMC);
1269 }
1270 else if(iS==2212){
1271 fContAllChargesMC->Fill(3,aodTrack->Pt(),v2mc,xMC);
587d006a 1272 }
1273 }
1274
afa8df58 1275 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1276
0f25ad32 1277 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 1278
587d006a 1279 Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1280 Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
a31d07c8 1281
243fbce7 1282 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
587d006a 1283 Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
afa8df58 1284 Float_t *probRead = fPID->GetProb();
1285 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1286 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2942f542 1287 Float_t x[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),1,static_cast<Float_t>(evPlAngV0[iV0]),static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to fill analysis v2 container
1288 Float_t x3[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),1,static_cast<Float_t>(evPlAngV0v3[iV0]),static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to fill analysis v3 container
0f25ad32 1289
1290 // in case fill DCA info
1291 if(fFillDCA){
1292 if(TMath::Abs(b[0]) > 0.1){
1293 x[5] = 1;
1294 x3[5] = 1;
1295 }
1296 if(TMath::Abs(b[0]) > 0.3){
1297 x[5] = 2;
1298 x3[5] = 2;
1299 }
1300 if(fIsMC && mcArray){
1301 if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
1302 x[5] += 3;
1303 x3[5] += 3;
1304 }
1305 }
1306 }
a640d925 1307
afa8df58 1308 // Fill no PID
587d006a 1309 if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
1310 if(fV3) contV0v3[iV0]->Fill(0,aodTrack->Pt(),v3V0,x3);
1311
afa8df58 1312
afa8df58 1313 Double_t dedxExp[8];
1314 Float_t tof = -1;
1315 Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1316 Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
78a615f6 1317
1318 Float_t nsigmaTPC[8];
1319 Float_t nsigmaTOF[8];
1320
afa8df58 1321 if(aodTrack->GetDetPid()){ // check the PID object is available
78a615f6 1322 for(Int_t iS=0;iS < 8;iS++){
afa8df58 1323 dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS);
78a615f6 1324 nsigmaTPC[iS] = (dedx - dedxExp[iS])/(dedxExp[iS]*0.07);
1325 // printf("TPC %i = %f (%f %f)\n",iS, nsigmaTPC[iS],dedx, dedxExp[iS]);
1326 }
587d006a 1327
afa8df58 1328 if(fPID->GetCurrentMask(1)){ // if TOF is present
1329 Float_t ptrack = aodTrack->P();
1330 tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack);
1331 aodTrack->GetIntegratedTimes(inttimes);
1332
1333 for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei
1334 inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]);
1335
78a615f6 1336 for(Int_t iS=0;iS<8;iS++){
1337 expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]);
1338 nsigmaTOF[iS] = (tof - inttimes[iS])/expTOFsigma[iS];
1339 // printf("TOF %i = %f\n",iS, nsigmaTOF[iS]);
1340 }
afa8df58 1341 }
1342 }
1343
1344 Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0];
1345 if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1346 else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1347 if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1348 else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1349
1350 Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0];
1351 if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1352 else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1353 if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1354 else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1355
1356 // variable to fill QA container
2942f542 1357 Float_t xQA[5] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt()), 0.0,deltaPhiV0,x[4]}; // v2
1358 Float_t xQA3[5] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt()), 0.0,deltaPhiV0v3,x[4]}; // v3
afa8df58 1359
0f25ad32 1360 // extra QA TProfiles
1361 if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
2942f542 1362 Float_t xQApid[2] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt())};
0f25ad32 1363 fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
1364 fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
1365 fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
1366 fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
1367 fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
1368 fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
1369 }
1370
78a615f6 1371 // QA fill
243fbce7 1372 if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
a640d925 1373 else{
243fbce7 1374 if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
1375 xQA[2] = prob[2];
1376 xQA3[2] = xQA[2];
a31d07c8 1377 if(fQAsw && fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
1378 if(fQAsw && fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
78a615f6 1379 }
243fbce7 1380 if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
1381 xQA[2] = prob[3];
1382 xQA3[2] = xQA[2];
a31d07c8 1383 if(fQAsw && fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
1384// if(fQAsw && fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);
78a615f6 1385 }
243fbce7 1386 if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
1387 xQA[2] = prob[4];
1388 xQA3[2] = xQA[2];
a31d07c8 1389 if(fQAsw && fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
1390// if(fQAsw && fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);
78a615f6 1391 }
243fbce7 1392 if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
1393 xQA[2] = prob[0];
1394 xQA3[2] = xQA[2];
a31d07c8 1395// if(fQAsw && fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
1396// if(fQAsw && fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);
78a615f6 1397 }
243fbce7 1398 if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
1399 xQA[2] = prob[5];
1400 xQA3[2] = xQA[2];
a31d07c8 1401 // if(fQAsw && fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
1402 // if(fQAsw && fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);
78a615f6 1403 }
243fbce7 1404 if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
1405 xQA[2] = prob[6];
1406 xQA3[2] = xQA[2];
a31d07c8 1407 // if(fQAsw && fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
1408 // if(fQAsw && fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);
78a615f6 1409 }
243fbce7 1410 if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
1411 xQA[2] = prob[7];
1412 xQA3[2] = xQA[2];
a31d07c8 1413 // if(fQAsw && fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
1414 // if(fQAsw && fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);
78a615f6 1415 }
78a615f6 1416 }
1417
afa8df58 1418 //pid selection
1419 if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!)
1420 else if(prob[2] > 0.6){ // pi
587d006a 1421 x[2] = prob[2];
587d006a 1422 x3[2] = x[2];
78a615f6 1423 if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
587d006a 1424 if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
a72e0394 1425 if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
0f25ad32 1426 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
1427 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1428 }
1429 }
1430 else if(prob[3] > 0.6){ // K
587d006a 1431 x[2] = prob[3];
587d006a 1432 x3[2] = x[2];
78a615f6 1433 if(TMath::Abs(nsigmaTPC[3]) < 5){
587d006a 1434 if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
a72e0394 1435 if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
0f25ad32 1436 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
1437 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1438 }
1439 }
1440 else if(prob[4] > 0.6){ // p
587d006a 1441 x[2] = prob[4];
587d006a 1442 x3[2] = x[2];
78a615f6 1443 if(TMath::Abs(nsigmaTPC[4]) < 5){
587d006a 1444 if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
a72e0394 1445 if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
0f25ad32 1446 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
1447 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1448 }
1449 }
1450 else if(prob[0] > 0.6){ // e
587d006a 1451 x[2] = prob[0];
587d006a 1452 x3[2] = x[2];
78a615f6 1453 if(TMath::Abs(nsigmaTPC[0]) < 5){
587d006a 1454 if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
a72e0394 1455 if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
0f25ad32 1456 if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
1457 else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1458 }
1459 }
587d006a 1460 else if(prob[1] > 0.6){ // mu
587d006a 1461 x[2] = prob[1];
587d006a 1462 x3[2] = x[2];
78a615f6 1463 if(TMath::Abs(nsigmaTPC[1]) < 5){
587d006a 1464 if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
a72e0394 1465 if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
587d006a 1466 }
1467 }
afa8df58 1468 else if(prob[5] > 0.6){ // d
587d006a 1469 x[2] = prob[5];
587d006a 1470 x3[2] = x[2];
78a615f6 1471 if(TMath::Abs(nsigmaTPC[5]) < 5){
587d006a 1472 if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
a72e0394 1473 if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
0f25ad32 1474 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
1475 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1476 }
1477 }
1478 else if(prob[6] > 0.6){ // t
587d006a 1479 x[2] = prob[6];
587d006a 1480 x3[2] = x[2];
78a615f6 1481 if(TMath::Abs(nsigmaTPC[6]) < 5){
587d006a 1482 if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
a72e0394 1483 if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
0f25ad32 1484 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
1485 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1486 }
1487 }
1488 else if(prob[7] > 0.6){ // He3
587d006a 1489 x[2] = prob[7];
587d006a 1490 x3[2] = x[2];
78a615f6 1491 if(TMath::Abs(nsigmaTPC[7]) < 5){
587d006a 1492 if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
a72e0394 1493 if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
0f25ad32 1494 if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
1495 else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
afa8df58 1496 }
afa8df58 1497 }
1498
587d006a 1499 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 1500 fPID->ResetDetOR(1); // exclude TOF from PID
1501 tofMismProb = 0;
1502
243fbce7 1503 fPID->ComputeProb(aodTrack,fOutputAOD);
587d006a 1504 dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
afa8df58 1505 probRead = fPID->GetProb();
1506
1507 fPID->SetDetOR(1); // include TOF for PID
1508 }
1509 Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities
1510
1511 //pid selection TPC S.A. with TOF matching
587d006a 1512 x[4]*=2; // set the mask to 2 id TOF is present
1513 if(x[4]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance
afa8df58 1514 else if(probTPC[2] > 0.6){ // pi
587d006a 1515 x[2] = probTPC[2];
1516 x3[2] = x[2];
78a615f6 1517 if(TMath::Abs(nsigmaTPC[2]) < 5){
587d006a 1518 if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
a72e0394 1519 if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
afa8df58 1520 }
1521 }
1522 else if(probTPC[3] > 0.6){ // K
587d006a 1523 x[2] = probTPC[3];
1524 x3[2] = x[2];
78a615f6 1525 if(TMath::Abs(nsigmaTPC[3]) < 5){
587d006a 1526 if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
a72e0394 1527 if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
afa8df58 1528 }
1529 }
1530 else if(probTPC[4] > 0.6){ // p
587d006a 1531 x[2] = probTPC[4];
1532 x3[2] = x[2];
78a615f6 1533 if(TMath::Abs(nsigmaTPC[4]) < 5){
587d006a 1534 if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
a72e0394 1535 if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
afa8df58 1536 }
1537 }
1538 else if(probTPC[0] > 0.6){ // e
587d006a 1539 x[2] = probTPC[0];
1540 x3[2] = x[2];
78a615f6 1541 if(TMath::Abs(nsigmaTPC[0]) < 5){
587d006a 1542 if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
a72e0394 1543 if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
587d006a 1544 }
1545 }
1546 else if(probTPC[1] > 0.6){ // mu
1547 x[2] = probTPC[1];
1548 x3[2] = x[2];
78a615f6 1549 if(TMath::Abs(nsigmaTPC[1]) < 5){
587d006a 1550 if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
a72e0394 1551 if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
afa8df58 1552 }
1553 }
1554 else if(probTPC[5] > 0.6){ // d
587d006a 1555 x[2] = probTPC[5];
1556 x3[2] = x[2];
78a615f6 1557 if(TMath::Abs(nsigmaTPC[5]) < 5){
587d006a 1558 if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
a72e0394 1559 if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
afa8df58 1560 }
1561 }
1562 else if(probTPC[6] > 0.6){ // t
587d006a 1563 x[2] = probTPC[6];
1564 x3[2] = x[2];
78a615f6 1565 if(TMath::Abs(nsigmaTPC[6]) < 5){
587d006a 1566 if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
a72e0394 1567 if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
afa8df58 1568 }
1569 }
1570 else if(probTPC[7] > 0.6){ // He3
587d006a 1571 x[2] = probTPC[7];
1572 x3[2] = x[2];
78a615f6 1573 if(TMath::Abs(nsigmaTPC[7]) < 5){
587d006a 1574 if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
a72e0394 1575 if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
afa8df58 1576 }
afa8df58 1577 }
1578 } // end side loop
1579 } // end track loop
1580
a31d07c8 1581 // my V0 loop
1582 for(Int_t imy=0;imy<fNK0s;imy++){
1583 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1584 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1585
1586 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1587 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1588
1589 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
2942f542 1590 Float_t x[6] = {static_cast<Float_t>(iC),-1/*my K0s are negative for convention*/,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1591 Float_t x3[6] = {static_cast<Float_t>(iC),-1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
a31d07c8 1592
1593 Float_t v2V0 = TMath::Cos(2*(fPhiK0s[imy] - evPlAngV0[iV0]));
1594 Float_t v3V0 = TMath::Cos(3*(fPhiK0s[imy] - evPlAngV0v3[iV0]));
1595 if(fV2) contV0[iV0]->Fill(9,fPtK0s[imy],v2V0,x);
1596 if(fV3) contV0v3[iV0]->Fill(9,fPtK0s[imy],v3V0,x3);
1597 }
1598 }
1599
a640d925 1600 // V0 loop
1601 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1602 AliAODv0 *myV0;
ae07aade 1603 Double_t /*dQT, dALPHA, dPT,*/ dMASS=0.0;
a640d925 1604 for (Int_t i=0; i!=nV0s; ++i) {
1605 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1606 if(!myV0) continue;
1607 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > fEtaCut) continue; // skipping low momentum
1608 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0);
1609 if(pass) {
1610 dMASS = myV0->MassK0Short();
1611 pass = 3;
a31d07c8 1612 fHK0sMass2->Fill(myV0->Pt(),dMASS);
a640d925 1613 }
a31d07c8 1614 if(TMath::Abs(dMASS-0.497)/0.005 > 3){
a640d925 1615 pass = PassesAODCuts(myV0,fOutputAOD,1);
1616 if(pass) dMASS = myV0->MassLambda();
1617 if(pass==2) dMASS = myV0->MassAntiLambda();
1618 }
a31d07c8 1619
a640d925 1620 if(pass){// 1 lambda, 2 antilambda, 3=K0s
ae07aade 1621 //dPT=myV0->Pt();
1622 //dQT=myV0->PtArmV0();
1623 //dALPHA=myV0->AlphaV0();
a640d925 1624
1625 Int_t iPos, iNeg;
1626 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0);
1627 if(iT->Charge()>0) {
1628 iPos = 0; iNeg = 1;
1629 } else {
1630 iPos = 1; iNeg = 0;
1631 }
a31d07c8 1632
1633 // check if one of the daugthers was already used
1634 if(pass == 3 && TMath::Abs(dMASS-0.497)/0.005 < 1){
1635 fHKsPhi->Fill(fZvtx, myV0->Phi());
1636 }
1637
1638 if(pass == 1000){ // disable
1639 Bool_t used = kFALSE;
1640 for(Int_t ii=0;ii<nusedForK0s;ii++){
1641 if(myV0->GetDaughter(iNeg) == usedForK0s1[ii] || myV0->GetDaughter(iPos) == usedForK0s2[ii]){
1642 used = kTRUE;
1643 }
1644 }
1645 if((!used) && nusedForK0s < 1000){
1646 nusedForK0s++;
1647 usedForK0s1[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iNeg);
1648 usedForK0s2[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iPos);
1649 printf("accepted\n");
1650 }
1651 else{
1652 dMASS = 0;
1653 printf("rejected\n");
1654 }
1655 }
1656
a640d925 1657 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1658 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1659
1660 // re-map the container in an array to do the analysis for V0A and V0C within a loop
1661 Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1662 AliFlowVZEROResults *contV0[2] = {fContAllChargesV0A,fContAllChargesV0C};
1663
1664 Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1665 AliFlowVZEROResults *contV0v3[2] = {fContAllChargesV0Av3,fContAllChargesV0Cv3};
1666
1667 for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1668
a31d07c8 1669 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)
1670
a640d925 1671 Float_t v2V0 = TMath::Cos(2*(myV0->Phi() - evPlAngV0[iV0]));
1672 Float_t v3V0 = TMath::Cos(3*(myV0->Phi() - evPlAngV0v3[iV0]));
1673
a31d07c8 1674 Float_t deltaphi = myV0->Phi()- evPlAngV0[iV0];
1675 if(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1676 if(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1677
2942f542 1678 Float_t x[6] = {static_cast<Float_t>(iC),1,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1679 Float_t x3[6] = {static_cast<Float_t>(iC),1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
a640d925 1680
1681 Float_t decaylength = myV0->DecayLengthXY(fOutputAOD->GetPrimaryVertex());
1682 // printf("decay length = %f\n",decaylength);
1683
1684 if(pass==2){ // anti-lambda charge = -1
1685 x[1] = -1;
1686 x3[1] = -1;
1687 }
1688
1689 if(decaylength < fMinDistV0) pass = 0;
1690 if(decaylength > fMaxDistV0) pass = 0;
1691
1692 Float_t nsigma = 0;
1693 if(pass < 3)
1694 nsigma = TMath::Abs(dMASS-1.116)/0.0016;
1695 else if(pass == 3)
1696 nsigma = TMath::Abs(dMASS-0.497)/0.005;
1697
1698 if(nsigma < 1)
1699 x[2] = 0.95;
1700 else if(nsigma < 2)
1701 x[2] = 0.85;
1702 else if(nsigma < 3)
1703 x[2] = 0.75;
1704 else if(nsigma < 4)
1705 x[2] = 0.65;
1706 else
1707 x[2] = 0.5;
1708
1709 x3[2] = x[2];
1710
1711 // Fill Container for lambda and Ks
a31d07c8 1712 if(fV2 && pass == 3 && x[2] > 0.6){
1713 contV0[iV0]->Fill(9,myV0->Pt(),v2V0,x);
1714 fHctauPtEP->Fill(myV0->Pt(),deltaphi,decaylength);//ciao
1715 if(myV0->Pt() < 1.1 && myV0->Pt() > 0.9) fHctauAt1EP->Fill(decaylength,deltaphi);
1716 }
a640d925 1717 if(fV3 && pass == 3 && x[2] > 0.6) contV0v3[iV0]->Fill(9,myV0->Pt(),v3V0,x3);
1718 if(fV2 && pass < 3 && x[2] > 0.6) contV0[iV0]->Fill(10,myV0->Pt(),v2V0,x);
1719 if(fV3 && pass < 3 && x[2] > 0.6) contV0v3[iV0]->Fill(10,myV0->Pt(),v3V0,x3);
1720
1721 if(pass < 3){ // lambda
1722 AliAODTrack* aodTrack = iT;
1723 if(pass==2) aodTrack=jT;
1724
1725 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1726 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1727
1728 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1729 Float_t *probRead = fPID->GetProb();
1730 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1731 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
a31d07c8 1732
1733 if(prob[4] < 0.61) prob[4] = 0.61;
a640d925 1734
2942f542 1735 Float_t xdec[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob[4],evPlAngV0[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to fill analysis v2 container
1736 Float_t xdec3[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob[4],evPlAngV0v3[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to fill analysis v3 container
a640d925 1737
1738 // Fill Container for (anti)proton from lambda
a31d07c8 1739 if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
a640d925 1740 if(fV2) contV0[iV0]->Fill(11,aodTrack->Pt(),v2V0,xdec);
1741 if(fV3) contV0v3[iV0]->Fill(11,aodTrack->Pt(),v3V0,xdec3);
1742 }
1743 }
a31d07c8 1744 else if(pass == 3){
1745 AliAODTrack* aodTrack = iT;
1746
1747 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1748 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1749
1750 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1751 Float_t *probRead = fPID->GetProb();
1752 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1753 Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1754
1755 if(prob[2] < 0.61) prob[2] = 0.61;
1756
2942f542 1757 Float_t xdec[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob[2],evPlAngV0[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to
1758 Float_t xdec3[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob[2],evPlAngV0v3[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb < 0.5),0}; // to
a31d07c8 1759
1760 if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1761 if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdec);
1762 if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdec3);
1763 }
1764
1765 aodTrack = jT;
1766 v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1767 v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1768
1769 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1770 Float_t *probRead2 = fPID->GetProb();
1771 Float_t prob2[8] = {probRead2[0],probRead2[1],probRead2[2],probRead2[3],probRead2[4],probRead2[5],probRead2[6],probRead2[7]};
1772 Float_t tofMismProb2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1773
1774 if(prob2[2] < 0.61) prob2[2] = 0.61;
1775
2942f542 1776 Float_t xdecB[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob2[2],evPlAngV0[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5),0}; // to
1777 Float_t xdecB3[6] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Charge()),prob2[2],evPlAngV0v3[iV0],static_cast<Float_t>(fPID->GetCurrentMask(1)&&tofMismProb2 < 0.5),0}; // to
a31d07c8 1778
1779 if(nsigma < 2 && xdecB[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1780 if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdecB);
1781 if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdecB3);
1782 }
1783 }
a640d925 1784 }
1785
1786 }
1787 } // end loop on V0
1788
1789
afa8df58 1790 // Fill EP distribution histograms
1791 if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
1792 if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
1793
1794 if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
1795 if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
1796
afa8df58 1797 // Fill histograms needed for resolution evaluation
1798 if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
1799 if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
1800 if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
1801
1802 if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
1803 if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
1804 if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
1805 }
1806
a31d07c8 1807
1808
1809 // clean track array
1810 for(Int_t i=0;i < nusedForK0s;i++){
1811 usedForK0s1[i] = NULL;
1812 usedForK0s2[i] = NULL;
1813 }
afa8df58 1814}
1815
1816//_____________________________________________________________________________
1817Float_t AliAnalysisTaskVnV0::GetVertex(AliAODEvent* aod) const
1818{
1819
1820 Float_t zvtx = -999;
1821
1822 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
1823 if (!vtxAOD)
1824 return zvtx;
1825 if(vtxAOD->GetNContributors()>0)
1826 zvtx = vtxAOD->GetZ();
1827
1828 return zvtx;
1829}
1830//_____________________________________________________________________________
1831void AliAnalysisTaskVnV0::Terminate(Option_t *)
1832{
1833 // Terminate loop
1834 Printf("Terminate()");
1835}
1836//_____________________________________________________________________________
1837void AliAnalysisTaskVnV0::OpenInfoCalbration(Int_t run){
1838 TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1839 TFile *foadb = TFile::Open(oadbfilename.Data());
1840
1841 if(!foadb){
1842 printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1843 return;
1844 }
1845
1846 AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1847 if(!cont){
1848 printf("OADB object hMultV0BefCorr is not available in the file\n");
1849 return;
1850 }
1851
1852 if(!(cont->GetObject(run))){
1853 printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1854 run = 137366;
1855 }
1856 fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1857
1858 TF1 *fpol0 = new TF1("fpol0","pol0");
1859 fMultV0->Fit(fpol0,"","",0,31);
1860 fV0Cpol = fpol0->GetParameter(0);
1861 fMultV0->Fit(fpol0,"","",32,64);
1862 fV0Apol = fpol0->GetParameter(0);
1863
1864 for(Int_t iside=0;iside<2;iside++){
1865 for(Int_t icoord=0;icoord<2;icoord++){
a640d925 1866 for(Int_t i=0;i < 9;i++){
afa8df58 1867 char namecont[100];
1868 if(iside==0 && icoord==0)
d6a1c304 1869 snprintf(namecont,100,"hQxc2_%i",i);
afa8df58 1870 else if(iside==1 && icoord==0)
d6a1c304 1871 snprintf(namecont,100,"hQxa2_%i",i);
afa8df58 1872 else if(iside==0 && icoord==1)
d6a1c304 1873 snprintf(namecont,100,"hQyc2_%i",i);
afa8df58 1874 else if(iside==1 && icoord==1)
d6a1c304 1875 snprintf(namecont,100,"hQya2_%i",i);
afa8df58 1876
1877 cont = (AliOADBContainer*) foadb->Get(namecont);
1878 if(!cont){
1879 printf("OADB object %s is not available in the file\n",namecont);
1880 return;
1881 }
1882
1883 if(!(cont->GetObject(run))){
1884 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1885 run = 137366;
1886 }
1887 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1888 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1889
1890 //for v3
1891 if(iside==0 && icoord==0)
d6a1c304 1892 snprintf(namecont,100,"hQxc3_%i",i);
afa8df58 1893 else if(iside==1 && icoord==0)
d6a1c304 1894 snprintf(namecont,100,"hQxa3_%i",i);
afa8df58 1895 else if(iside==0 && icoord==1)
d6a1c304 1896 snprintf(namecont,100,"hQyc3_%i",i);
afa8df58 1897 else if(iside==1 && icoord==1)
d6a1c304 1898 snprintf(namecont,100,"hQya3_%i",i);
afa8df58 1899
1900 cont = (AliOADBContainer*) foadb->Get(namecont);
1901 if(!cont){
1902 printf("OADB object %s is not available in the file\n",namecont);
1903 return;
1904 }
1905
1906 if(!(cont->GetObject(run))){
1907 printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1908 run = 137366;
1909 }
1910 fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1911 fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1912
1913 }
1914 }
1915 }
1916}
a640d925 1917//=======================================================================
1918Int_t AliAnalysisTaskVnV0::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1919{
a31d07c8 1920 Int_t set = 0;
a640d925 1921 Float_t fV0Cuts[9];
1922 // defines cuts to be used
1923 // fV0Cuts[9] dl dca ctp d0 d0d0 qt minEta maxEta PID
1924 switch(set) {
1925 case(0): // No cuts
1926 fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1927 fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1928 fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 0;
1929 break;
1930 case(1): // Tight cuts
1931 fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1932 fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1933 fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 0;
1934 break;
1935 case(2): // Tight cuts + PID
1936 fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1937 fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1938 fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 1;
1939 break;
1940 case(3): // No cuts + PID
1941 fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1942 fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1943 fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 1;
1944 break;
1945 }
1946
1947 // daughter cuts
1948 if(! fCutsDaughter){
1949 fCutsDaughter = new AliESDtrackCuts(Form("daughter_cuts_%s","ESD") );
1950 fCutsDaughter->SetPtRange(0.2,10.0);
1951 fCutsDaughter->SetEtaRange(-0.8, 0.8 );
1952 fCutsDaughter->SetMinNClustersTPC(80);
1953 fCutsDaughter->SetMaxChi2PerClusterTPC(4.0);
1954 fCutsDaughter->SetRequireTPCRefit(kTRUE);
1955 fCutsDaughter->SetAcceptKinkDaughters(kFALSE);
1956 }
1957
1958 if (myV0->GetOnFlyStatus() ) return 0;
1959 //the following is needed in order to evualuate track-quality
1960 AliAODTrack *iT, *jT;
1961 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1962 Double_t pos[3],cov[6];
1963 vV0s->GetXYZ(pos);
1964 vV0s->GetCovarianceMatrix(cov);
1965 const AliESDVertex vESD(pos,cov,100.,100);
1966 // TESTING CHARGE
1967 int iPos, iNeg;
1968 iT=(AliAODTrack*) myV0->GetDaughter(0);
1969 if(iT->Charge()>0) {
1970 iPos = 0; iNeg = 1;
1971 } else {
1972 iPos = 1; iNeg = 0;
1973 }
1974 // END OF TEST
1975
1976 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1977 AliESDtrack ieT( iT );
1978 ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1979 ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1980 ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1981 ieT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1982 if (!fCutsDaughter->IsSelected( &ieT ) ) return 0;
1983
1984 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1985 AliESDtrack jeT( jT );
1986 jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1987 jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1988 jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1989 jeT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1990 if (!fCutsDaughter->IsSelected( &jeT ) ) return 0;
1991
1992 Double_t pvertex[3];
1993 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1994 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1995 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1996 Double_t dDL=myV0->DecayLengthV0( pvertex );
1997 Double_t dDCA=myV0->DcaV0Daughters();
1998 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1999 Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2000 Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2001 Double_t dD0D0=dD0P*dD0M;
2002 Double_t dQT=myV0->PtArmV0();
2003 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
2004 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
2005// Double_t dPT=myV0->Pt();
2006 Double_t dETA=myV0->Eta();
2007 Int_t passes = 1;
2008 if(dDL <fV0Cuts[0]) passes = 0;
2009 if(dDCA >fV0Cuts[1]) passes = 0;
2010 if(dCTP <fV0Cuts[2]) passes = 0;
2011 if(TMath::Abs(dD0P) <fV0Cuts[3]) passes = 0;
2012 if(TMath::Abs(dD0M) <fV0Cuts[3]) passes = 0;
2013 if(dD0D0>fV0Cuts[4]) passes = 0;
2014 if(dETA <fV0Cuts[6]) passes = 0;
2015 if(dETA >fV0Cuts[7]) passes = 0;
2016 if(specie==0) if(dQT<fV0Cuts[5]) passes = 0;
2017 if(specie==1&&passes==1&&dALPHA<0) passes = 2; // antilambda
a31d07c8 2018
2019
2020// if(jT->Pt() < 0.5*myV0->Pt() || iT->Pt() < 0.5*myV0->Pt()) passes = 0;
2021
2022
2023 // additional cut
2024// if(!(iT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2025// if(!(jT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2026
2027// if(!(iT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2028// if(!(jT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2029
2030// if(!(iT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2031// if(!(jT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2032
2033 Bool_t trkFlag = iT->TestFilterBit(1); // TPC only tracks (4,global track)
2034 Bool_t trkFlag2 = jT->TestFilterBit(1); // TPC only tracks (4,global track)
2035
2036 if(!trkFlag) passes = 0;
2037 if(!trkFlag2) passes = 0;
2038
a640d925 2039 if(passes&&fV0Cuts[8]) {
2040
2041 Double_t dedxExp[8];
2042 fPID->ComputeProb(iT,tAOD); // compute Bayesian probabilities
2043 Float_t nsigmaTPC[8];
a31d07c8 2044
ae07aade 2045// Int_t tofMatch=0;
2046// Int_t tofMatch2=0;
a31d07c8 2047
a640d925 2048 if(iT->GetDetPid()){ // check the PID object is available
2049 for(Int_t iS=0;iS < 8;iS++){
2050 dedxExp[iS] = fPID->GetExpDeDx(iT,iS);
2051 nsigmaTPC[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2052 }
2053 }
2054 else{
2055 for(Int_t iS=0;iS < 8;iS++)
2056 nsigmaTPC[iS] = 10;
2057 }
2058
a31d07c8 2059
ae07aade 2060// if(fPID->GetCurrentMask(1)) // if TOF is present
2061// tofMatch = 1;
a31d07c8 2062
2063// Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2064
2065 Float_t *probRead = fPID->GetProb();
2066 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2067
a640d925 2068 fPID->ComputeProb(jT,tAOD); // compute Bayesian probabilities
2069 Float_t nsigmaTPC2[8];
2070 if(jT->GetDetPid()){ // check the PID object is available
2071 for(Int_t iS=0;iS < 8;iS++){
2072 dedxExp[iS] = fPID->GetExpDeDx(jT,iS);
2073 nsigmaTPC2[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2074 }
2075 }
2076 else{
2077 for(Int_t iS=0;iS < 8;iS++)
2078 nsigmaTPC2[iS] = 10;
2079 }
2080
ae07aade 2081// if(fPID->GetCurrentMask(1)) // if TOF is present
2082// tofMatch2 = 1;
a31d07c8 2083
2084// Float_t tofMismProbMC2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2085
2086 probRead = fPID->GetProb();
2087 Float_t prob2[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2088
a640d925 2089 if(jT->GetTPCNcls() < fNcluster) passes = 0;
2090 else if(iT->GetTPCNcls() < fNcluster) passes = 0;
2091
a31d07c8 2092// if(! (tofMatch && tofMatch2)) passes = 0;
2093
2094 /*
2095 Float_t dMASS = myV0->MassK0Short();
2096 Float_t nsigmaMass = TMath::Abs(dMASS-0.497)/0.005;
2097 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);
2098 */
2099
a640d925 2100 switch(specie) {
2101 case 0: // K0 PID
a31d07c8 2102 if(0){
2103 if( ((jT->GetTPCmomentum()<15) &&
2104 (TMath::Abs(nsigmaTPC2[2])>3.)) || prob2[2] < 0.9)
2105 passes = 0;
2106 if( ((iT->GetTPCmomentum()<15) &&
2107 (TMath::Abs(nsigmaTPC[2])>3.))|| prob[2] < 0.9 )
2108 passes = 0;
2109 }
a640d925 2110 break;
2111 case 1: // Lambda PID i==pos j ==neg
2112 if(passes==1) {
2113 if( (iT->GetTPCmomentum()<15) &&
2114 (TMath::Abs(nsigmaTPC[4])>3.) )
2115 passes = 0;
2116 if( (jT->GetTPCmomentum()<15) &&
2117 (TMath::Abs(nsigmaTPC2[2])>3.) )
2118 passes = 0;
2119 }
2120 if(passes==2) {
2121 if( (iT->GetTPCmomentum()<15) &&
2122 (TMath::Abs(nsigmaTPC[2])>3.) )
2123 passes = 0;
2124 if( (jT->GetTPCmomentum()<15) &&
2125 (TMath::Abs(nsigmaTPC2[4])>3.) )
2126 passes = 0;
2127 }
2128 break;
2129 }
2130 }
2131 return passes;
2132}
a31d07c8 2133//=======================================================================
2134void AliAnalysisTaskVnV0::SelectK0s(){
2135 fNK0s=0;
2136 fNpiPos=0;
2137 fNpiNeg=0;
2138
2139 if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAng2,1.0); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
2140
2141 // fill pion stacks
2142 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
2143 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
2144 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
2145
2146 if (!aodTrack){
a31d07c8 2147 continue;
2148 }
2149
2150 Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
2151// trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
2152
2153 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
2154 continue;
2155 }
2156
2157 Double_t b[2] = {-99., -99.};
2158 Double_t bCov[3] = {-99., -99., -99.};
d8884981 2159
4d0ef7e4 2160 AliAODTrack param(*aodTrack);
2161 if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
a31d07c8 2162 continue;
d8884981 2163 }
4d0ef7e4 2164
a31d07c8 2165 if(TMath::Abs(b[0]) < 0.5/aodTrack->Pt()) continue;
2166
2167 fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
2168 Float_t *probRead = fPID->GetProb();
2169 Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2170 // Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2171
2172 Int_t charge = aodTrack->Charge();
2173 if(prob[2] > 0.9){
2174 if(charge > 0){
2175 fIPiPos[fNpiPos] = iT;
2176 fNpiPos++;
2177 }
2178 else{
2179 fIPiNeg[fNpiNeg] = iT;
2180 fNpiNeg++;
2181 }
2182 }
2183 }
2184
2185 for(Int_t i=0;i < fNpiPos;i++){
2186 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
2187 AliESDtrack pipE(pip);
2188
2189 for(Int_t j=0;j < fNpiNeg;j++){
2190 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
2191 AliESDtrack pinE(pin);
2192
2193 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
2194
2195 Double_t pPos[3];
2196 Double_t pNeg[3];
2197 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
2198 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
2199
2200 Float_t length = (xp+xn)*0.5;
2201
2202 Float_t pxs = pPos[0] + pNeg[0];
2203 Float_t pys = pPos[1] + pNeg[1];
2204 Float_t pzs = pPos[2] + pNeg[2];
2205 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);
2206
2207 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
2208 Float_t phi = TMath::ATan2(pys,pxs);
2209 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
2210
2211 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
2212
2213 if(mindist < 0.2&& length > 1 && length < 25){
2214 fHK0sMass->Fill(pt,mass);
2215
2216 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);
2217 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);
2218
2219 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
2220 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
2221
2222 fHK0vsLambda->Fill(mass,TMath::Min(massaL,massaAL));
2223
2224 if(TMath::Abs(mass-0.497)/0.005 < 1 && massaL > 1.15 && massaAL > 1.15){
2225 fPhiK0s[fNK0s] = phi;
2226 fPtK0s[fNK0s] = pt;
2227 fNK0s++;
2228 }
2229 }
2230 }
2231 }
2232}