AliPhysics  b97afa6 (b97afa6)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliAnalysisTaskVnV0.cxx
Go to the documentation of this file.
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"
18 #include "AliGenHijingEventHeader.h"
19 #include "AliMCEvent.h"
20 #include "AliAODMCHeader.h"
21 #include "AliAODMCParticle.h"
22 #include "TChain.h"
23 #include "AliESDtrackCuts.h"
24 #include "AliESDVertex.h"
25 #include "AliEventplane.h"
26 #include "TProfile2D.h"
27 
28 // STL includes
29 //#include <iostream>
30 //using namespace std;
31 
32 ClassImp(AliAnalysisTaskVnV0)
33 Bool_t AliAnalysisTaskVnV0::fgIsPsiComputed = kFALSE;
34 Float_t AliAnalysisTaskVnV0::fgPsi2v0a=999.;
35 Float_t AliAnalysisTaskVnV0::fgPsi2v0c=999.;
36 Float_t AliAnalysisTaskVnV0::fgPsi2tpc=999.;
37 Float_t AliAnalysisTaskVnV0::fgPsi3v0a=999.;
38 Float_t AliAnalysisTaskVnV0::fgPsi3v0c=999.;
39 Float_t AliAnalysisTaskVnV0::fgPsi3tpc=999.;
40 Float_t AliAnalysisTaskVnV0::fgPsi2v0aMC=999.;
41 Float_t AliAnalysisTaskVnV0::fgPsi2v0cMC=999.;
42 Float_t AliAnalysisTaskVnV0::fgPsi2tpcMC=999.;
43 Float_t AliAnalysisTaskVnV0::fgPsi3v0aMC=999.;
44 Float_t AliAnalysisTaskVnV0::fgPsi3v0cMC=999.;
45 Float_t AliAnalysisTaskVnV0::fgPsi3tpcMC=999.;
46 
47 //_____________________________________________________________________________
50  fVtxCut(10.0), // cut on |vertex| < fVtxCut
51  fEtaCut(0.8), // cut on |eta| < fEtaCut
52  fMinPt(0.15), // cut on pt > fMinPt
53  fMinDistV0(0),
54  fMaxDistV0(100),
55  fV2(kTRUE),
56  fV3(kTRUE),
57  fIsMC(kFALSE),
58  fQAsw(kFALSE),
59  fIsAfter2011(kFALSE),
60  fRun(-1),
61  fNcluster(70),
62  fList(new TList()),
63  fList2(new TList()),
64  fList3(new TList()),
65  fList4(new TList()),
66  fMultV0(NULL),
67  fV0Cpol(100),
68  fV0Apol(100),
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),
79  fQA(NULL),
80  fQA2(NULL),
81  fQAv3(NULL),
82  fQA2v3(NULL),
83  fPID(new AliFlowBayesianPID()),
84  fTree(NULL),
85  fCentrality(-1),
86  evPlAngV0ACor2(0),
87  evPlAngV0CCor2(0),
88  evPlAng2(0),
89  evPlAngV0ACor3(0),
90  evPlAngV0CCor3(0),
91  evPlAng3(0),
92  fContAllChargesV0A(NULL),
93  fContAllChargesV0C(NULL),
94  fContAllChargesV0Av3(NULL),
95  fContAllChargesV0Cv3(NULL),
96  fContAllChargesMC(NULL),
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),
106  fContAllChargesMCCv3(NULL),
107  fFillDCA(kFALSE),
108  fContQApid(NULL),
109  fModulationDEDx(kFALSE),
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),
121  fCutsDaughter(NULL)
122 {
123  // Default constructor (should not be used)
124  fList->SetName("resultsV2");
125  fList2->SetName("resultsV3");
126  fList3->SetName("resultsMC");
127  fList4->SetName("QA");
128 
129  fList->SetOwner(kTRUE);
130  fList2->SetOwner(kTRUE);
131  fList3->SetOwner(kTRUE);
132  fList4->SetOwner(kTRUE);
133 
134  fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
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  }
142 }
143 
144 //______________________________________________________________________________
146  AliAnalysisTaskSE(name),
147  fVtxCut(10.0), // cut on |vertex| < fVtxCut
148  fEtaCut(0.8), // cut on |eta| < fEtaCut
149  fMinPt(0.15), // cut on pt > fMinPt
150  fMinDistV0(0),
151  fMaxDistV0(100),
152  fV2(kTRUE),
153  fV3(kTRUE),
154  fIsMC(kFALSE),
155  fQAsw(kFALSE),
156  fIsAfter2011(kFALSE),
157  fRun(-1),
158  fNcluster(70),
159  fList(new TList()),
160  fList2(new TList()),
161  fList3(new TList()),
162  fList4(new TList()),
163  fMultV0(NULL),
164  fV0Cpol(100),
165  fV0Apol(100),
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),
176  fQA(NULL),
177  fQA2(NULL),
178  fQAv3(NULL),
179  fQA2v3(NULL),
180  fPID(new AliFlowBayesianPID()),
181  fTree(NULL),
182  fCentrality(-1),
183  evPlAngV0ACor2(0),
184  evPlAngV0CCor2(0),
185  evPlAng2(0),
186  evPlAngV0ACor3(0),
187  evPlAngV0CCor3(0),
188  evPlAng3(0),
189  fContAllChargesV0A(NULL),
190  fContAllChargesV0C(NULL),
191  fContAllChargesV0Av3(NULL),
192  fContAllChargesV0Cv3(NULL),
193  fContAllChargesMC(NULL),
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),
203  fContAllChargesMCCv3(NULL),
204  fFillDCA(kFALSE),
205  fContQApid(NULL),
206  fModulationDEDx(kFALSE),
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),
218  fCutsDaughter(NULL)
219 {
220 
221  DefineOutput(1, TList::Class());
222  DefineOutput(2, TList::Class());
223  DefineOutput(3, TList::Class());
224  DefineOutput(4, TList::Class());
225 
226  // Output slot #1 writes into a TTree
227  fList->SetName("resultsV2");
228  fList2->SetName("resultsV3");
229  fList3->SetName("resultsMC");
230  fList4->SetName("QA");
231 
232  fList->SetOwner(kTRUE);
233  fList2->SetOwner(kTRUE);
234  fList3->SetOwner(kTRUE);
235  fList4->SetOwner(kTRUE);
236 
237  fPID->SetNewTrackParam(); // Better tuning for TOF PID tracking effect in LHC10h
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  }
245 }
246 //_____________________________________________________________________________
248 {
249 
250 }
251 
252 //______________________________________________________________________________
254 {
255 
256  if(fIsMC) fPID->SetMC(kTRUE);
257 
258 
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};
273  const Int_t nCentrTOF = nCentrBin;
274  const Int_t nPsiTOF = 10;
275  const Int_t nChargeBinsTOFres = 2;
276  const Int_t nCentrTOFres = nCentrBin;
277  const Int_t nProbTOFres = 4;
278  const Int_t nPsiTOFres = 10;
279  const Int_t nMaskPID = 3;
280 
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};
292  Int_t binsTOFmc[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,2};
293  Int_t binsTOFmcPureMC[5] = {nCentrTOFres,nChargeBinsTOFres,1,nPsiTOFres,1};
294 
295  // v2 container
296  fContAllChargesV0A = new AliFlowVZEROResults("v2A",6,binsTOF);
297  fContAllChargesV0A->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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
302  fContAllChargesV0A->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
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");
308  fContAllChargesV0A->SetVarName(5,"DCAbin");
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);
318  if(fV2) fContAllChargesV0A->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
319  if(fV2) fContAllChargesV0A->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
320  if(fV2) fContAllChargesV0A->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
321  if(fV2) fContAllChargesV0A->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
322 
323  fContAllChargesV0C = new AliFlowVZEROResults("v2C",6,binsTOF);
324  fContAllChargesV0C->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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
329  fContAllChargesV0C->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
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");
335  fContAllChargesV0C->SetVarName(5,"DCAbin");
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);
345  if(fV2) fContAllChargesV0C->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
346  if(fV2) fContAllChargesV0C->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
347  if(fV2) fContAllChargesV0C->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
348  if(fV2) fContAllChargesV0C->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
349 
352 
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 
357  if(fIsMC && fV2){
358  fContAllChargesMC = new AliFlowVZEROResults("v2mc",5,binsTOFmc);
359  fContAllChargesMC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
376 
377  fContAllChargesMCA = new AliFlowVZEROResults("v2mcA",5,binsTOFmcPureMC);
378  fContAllChargesMCA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
397  fContAllChargesMCC->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
414  }
415 
416  // v3 container
417  fContAllChargesV0Av3 = new AliFlowVZEROResults("v3A",6,binsTOF);
418  fContAllChargesV0Av3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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
422  fContAllChargesV0Av3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
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");
429  fContAllChargesV0Av3->SetVarName(5,"DCAbin");
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);
439  if(fV3) fContAllChargesV0Av3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
440  if(fV3) fContAllChargesV0Av3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
441  if(fV3) fContAllChargesV0Av3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
442  if(fV3) fContAllChargesV0Av3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
443 
444  fContAllChargesV0Cv3 = new AliFlowVZEROResults("v3C",6,binsTOF);
445  fContAllChargesV0Cv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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
450  fContAllChargesV0Cv3->SetVarRange(5,-0.5,nDCABin-0.5); // DCA mask
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");
456  fContAllChargesV0Cv3->SetVarName(5,"DCAbin");
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);
466  if(fV3) fContAllChargesV0Cv3->AddSpecies("Ks",nPtBinsTOF,binsPtTOF);
467  if(fV3) fContAllChargesV0Cv3->AddSpecies("Lambda",nPtBinsTOF,binsPtTOF);
468  if(fV3) fContAllChargesV0Cv3->AddSpecies("pFromLambda",nPtBinsTOF,binsPtTOF);
469  if(fV3) fContAllChargesV0Cv3->AddSpecies("piFromK",nPtBinsTOF,binsPtTOF);
470 
473 
474  if(fIsMC && fV3){
475  fContAllChargesMCAv3 = new AliFlowVZEROResults("v3mcA",5,binsTOFmcPureMC);
476  fContAllChargesMCAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
493 
494  fContAllChargesMCCv3 = new AliFlowVZEROResults("v3mcC",5,binsTOFmcPureMC);
495  fContAllChargesMCCv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
512  }
513 
514  // TProfile for resolutions 3 subevents (V0A, V0C, TPC)
515  // v2
516  fHResTPCv0A2 = new TProfile("hResTPCv0A2","",nCentrBin,0,nCentrBin);
517  fHResTPCv0C2 = new TProfile("hResTPCv0C2","",nCentrBin,0,nCentrBin);
518  fHResv0Cv0A2 = new TProfile("hResv0Cv0A2","",nCentrBin,0,nCentrBin);
519 
520  fList->Add(fHResTPCv0A2);
521  fList->Add(fHResTPCv0C2);
522  fList->Add(fHResv0Cv0A2);
523 
524  // v3
525  fHResTPCv0A3 = new TProfile("hResTPCv0A3","",nCentrBin,0,nCentrBin);
526  fHResTPCv0C3 = new TProfile("hResTPCv0C3","",nCentrBin,0,nCentrBin);
527  fHResv0Cv0A3 = new TProfile("hResv0Cv0A3","",nCentrBin,0,nCentrBin);
528 
529  fList2->Add(fHResTPCv0A3);
530  fList2->Add(fHResTPCv0C3);
531  fList2->Add(fHResv0Cv0A3);
532 
533  // MC as in the dataEP resolution (but using MC tracks)
534  if(fIsMC && fV3){
535  fHResMA2 = new TProfile("hResMA2","",nCentrBin,0,nCentrBin);
536  fHResMC2 = new TProfile("hResMC2","",nCentrBin,0,nCentrBin);
537  fHResAC2 = new TProfile("hResAC2","",nCentrBin,0,nCentrBin);
538  fList3->Add(fHResMA2);
539  fList3->Add(fHResMC2);
540  fList3->Add(fHResAC2);
541  }
542  if(fIsMC && fV3){
543  fHResMA3 = new TProfile("hResMA3","",nCentrBin,0,nCentrBin);
544  fHResMC3 = new TProfile("hResMC3","",nCentrBin,0,nCentrBin);
545  fHResAC3 = new TProfile("hResAC3","",nCentrBin,0,nCentrBin);
546  fList3->Add(fHResMA3);
547  fList3->Add(fHResMC3);
548  fList3->Add(fHResAC3);
549  }
550 
551 
552  // V0A and V0C event plane distributions
553  //v2
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);
556 
557  //v3
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);
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  }
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 // }
573 
574  const Int_t nDeltaPhi = 5;
575  const Int_t nDeltaPhiV3 = 7;
576 
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);
582  fQA->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
601  fQA2->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
620  fQAv3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
639  fQA2v3->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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);
656 
657  fList->Add(fPhiRPv0A);
658  fList->Add(fPhiRPv0C);
659 
660  if(fQAsw && fV2){
661  fList4->Add(fQA);
662  fList4->Add(fQA2);
663  }
664 
665  fList2->Add(fPhiRPv0Av3);
666  fList2->Add(fPhiRPv0Cv3);
667 
668  if(fQAsw && fV3){
669  fList4->Add(fQAv3);
670  fList4->Add(fQA2v3);
671  }
672 
673  // fList->Add(fTree); // comment if not needed
674 
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);
738  fContQApid->SetVarRange(0,-0.5,nCentrBin-0.5); // centrality
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 
750  printf("Output creation ok!!\n\n\n\n");
751 
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 
768  // Post output data.
769  if(fV2) PostData(1, fList);
770  if(fV3) PostData(2, fList2);
771  if(fIsMC) PostData(3, fList3);
772  if(fQAsw) PostData(4, fList4);
773 
774 }
775 
776 //______________________________________________________________________________
778 {
779  // Main loop
780  // Called for each event
781 
782  fgIsPsiComputed = kFALSE;
783  fgPsi2v0a=999.;
784  fgPsi2v0c=999.;
785  fgPsi2tpc=999.;
786  fgPsi3v0a=999.;
787  fgPsi3v0c=999.;
788  fgPsi3tpc=999.;
789  fgPsi2v0aMC=999.;
790  fgPsi2v0cMC=999.;
791  fgPsi2tpcMC=999.;
792  fgPsi3v0aMC=999.;
793  fgPsi3v0cMC=999.;
794  fgPsi3tpcMC=999.;
795 
796  fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
797  if(!fOutputAOD){
798  Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
799  this->Dump();
800  return;
801  }
802 
803  Int_t run = fOutputAOD->GetRunNumber();
804 
805  if(run != fRun){
806  // Load the calibrations run dependent
807  if(! fIsAfter2011) OpenInfoCalbration(run);
808  fRun=run;
809  }
810 
811  fZvtx = GetVertex(fOutputAOD);
812 
813 
814 
815  //Get the MC object
816  if(fIsMC){
817  AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
818  if (!mcHeader) {
819  AliError("Could not find MC Header in AOD");
820  return;
821  }
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 
844  if (TMath::Abs(fZvtx) < fVtxCut) {
845  //Centrality
846  Float_t v0Centr = -10.;
847  Float_t trkCentr = -10.;
848  AliCentrality *centrality = fOutputAOD->GetCentrality();
849  if (centrality){
850 // printf("v0centr = %f -- tpccnetr%f\n",centrality->GetCentralityPercentile("V0M"),centrality->GetCentralityPercentile("TRK"));
851  v0Centr = centrality->GetCentralityPercentile("V0M");
852  trkCentr = centrality->GetCentralityPercentile("TRK");
853  //changed
854  }
855 
856  if(TMath::Abs(v0Centr - trkCentr) < 5.0 && v0Centr > 0){ // consistency cut on centrality selection
857  fPID->SetDetResponse(fOutputAOD, v0Centr); // Set the PID object for each event!!!!
858  Analyze(fOutputAOD,v0Centr); // Do analysis!!!!
859 
860  fCentrality = v0Centr;
861  if(fV2) fTree->Fill();
862  }
863  }
864 
865 }
866 
867 //________________________________________________________________________
869 {
870 
871  Int_t nusedForK0s=0;
872  AliAODTrack *usedForK0s1[1000];
873  AliAODTrack *usedForK0s2[1000];
874 
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;
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
885  fgIsPsiComputed = kTRUE;
886 
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;
897 
898  Int_t iCcal = iC;
899 
900 /*
901  if(nCentrBin==16){
902  iC = v0Centr/5;
903  if(iC >= nCentrBin) iC = nCentrBin-1;
904  }
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 
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;
933 
934  Int_t nAODTracks = aodEvent->GetNumberOfTracks();
935 
936  AliAODMCHeader *mcHeader = NULL;
937  TClonesArray *mcArray = NULL;
938  Float_t evplaneMC = 0;
939  if(fIsMC){
940  mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
941 
942  if (mcHeader) {
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();
946  mcArray = (TClonesArray*)fOutputAOD->GetList()->FindObject(AliAODMCParticle::StdBranchName());
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){
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.;
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])));
986  fgPsi2v0aMC = EvPlaneMCV2[0];
987  fgPsi2v0cMC = EvPlaneMCV2[1];
988  fgPsi2tpcMC = EvPlaneMCV2[2];
989  }
990  if(fV3){
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.;
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])));
997  fgPsi3v0aMC = EvPlaneMCV3[0];
998  fgPsi3v0cMC = EvPlaneMCV3[1];
999  fgPsi3tpcMC = EvPlaneMCV3[2];
1000  }
1001 
1002  // flow A and C side
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};
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 
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);
1033 
1034  if(iS==11){
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);
1039  }
1040  else if(iS==13){
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);
1045  }
1046  else if(iS==211){
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);
1051  }
1052  else if(iS==321){
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);
1057  }
1058  else if(iS==2212){
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);
1063  }
1064  }
1065  }
1066  }
1067  }
1068 
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 = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(iT));
1076  if(!aodTrack) AliFatal("Not a standard AOD");
1077 
1078  if (!aodTrack){
1079  continue;
1080  }
1081 
1082  Bool_t trkFlag = aodTrack->TestFilterBit(1);
1083 
1084  if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
1085  continue;
1086 
1087  Double_t b[2] = {-99., -99.};
1088  Double_t bCov[3] = {-99., -99., -99.};
1089 
1090 
1091  AliAODTrack param(*aodTrack);
1092  if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
1093  continue;
1094  }
1095 
1096  if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
1097  continue;
1098 
1099  Qx2 += TMath::Cos(2*aodTrack->Phi());
1100  Qy2 += TMath::Sin(2*aodTrack->Phi());
1101  Qx3 += TMath::Cos(3*aodTrack->Phi());
1102  Qy3 += TMath::Sin(3*aodTrack->Phi());
1103 
1104  }
1105 
1106  evPlAng2 = TMath::ATan2(Qy2, Qx2)/2.;
1107  evPlAng3 = TMath::ATan2(Qy3, Qx3)/3.;
1108 
1109  fgPsi2tpc = evPlAng2;
1110  fgPsi3tpc = evPlAng3;
1111 
1112  SelectK0s();
1113 
1114  //V0 info
1115  AliAODVZERO* aodV0 = aodEvent->GetVZEROData();
1116 
1117  for (Int_t iv0 = 0; iv0 < 64; iv0++) {
1118  Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
1119  Float_t multv0 = aodV0->GetMultiplicity(iv0);
1120 
1121  if(! fIsAfter2011){
1122  if(! fIsMC){
1123  if (iv0 < 32){ // V0C
1124  Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1125  Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1126  Qxc3 += TMath::Cos(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1127  Qyc3 += TMath::Sin(3*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1128  } else { // V0A
1129  Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1130  Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1131  Qxa3 += TMath::Cos(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1132  Qya3 += TMath::Sin(3*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
1133  }
1134  }
1135  else{
1136  if (iv0 < 32){ // V0C
1137  Qxc2 += TMath::Cos(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1138  Qyc2 += TMath::Sin(2*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1139  Qxc3 += TMath::Cos(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1140  Qyc3 += TMath::Sin(3*phiV0) * multv0;//*fV0Cpol/fMultV0->GetBinContent(iv0+1);
1141  } else { // V0A
1142  Qxa2 += TMath::Cos(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1143  Qya2 += TMath::Sin(2*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1144  Qxa3 += TMath::Cos(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1145  Qya3 += TMath::Sin(3*phiV0) * multv0;//*fV0Apol/fMultV0->GetBinContent(iv0+1);
1146  }
1147  }
1148  }
1149  }
1150 
1151  //grab for each centrality the proper histo with the Qx and Qy to do the recentering
1152  Double_t Qxamean2 = fMeanQ[iCcal][1][0];
1153  Double_t Qxarms2 = fWidthQ[iCcal][1][0];
1154  Double_t Qyamean2 = fMeanQ[iCcal][1][1];
1155  Double_t Qyarms2 = fWidthQ[iCcal][1][1];
1156  Double_t Qxamean3 = fMeanQv3[iCcal][1][0];
1157  Double_t Qxarms3 = fWidthQv3[iCcal][1][0];
1158  Double_t Qyamean3 = fMeanQv3[iCcal][1][1];
1159  Double_t Qyarms3 = fWidthQv3[iCcal][1][1];
1160 
1161  Double_t Qxcmean2 = fMeanQ[iCcal][0][0];
1162  Double_t Qxcrms2 = fWidthQ[iCcal][0][0];
1163  Double_t Qycmean2 = fMeanQ[iCcal][0][1];
1164  Double_t Qycrms2 = fWidthQ[iCcal][0][1];
1165  Double_t Qxcmean3 = fMeanQv3[iCcal][0][0];
1166  Double_t Qxcrms3 = fWidthQv3[iCcal][0][0];
1167  Double_t Qycmean3 = fMeanQv3[iCcal][0][1];
1168  Double_t Qycrms3 = fWidthQv3[iCcal][0][1];
1169 
1170  Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
1171  Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
1172  Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
1173  Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
1174  Double_t QxaCor3 = (Qxa3 - Qxamean3)/Qxarms3;
1175  Double_t QyaCor3 = (Qya3 - Qyamean3)/Qyarms3;
1176  Double_t QxcCor3 = (Qxc3 - Qxcmean3)/Qxcrms3;
1177  Double_t QycCor3 = (Qyc3 - Qycmean3)/Qycrms3;
1178 
1179  if(! fIsAfter2011){
1180  if(! fIsMC){
1181  evPlAngV0ACor2 = TMath::ATan2(QyaCor2, QxaCor2)/2.;
1182  evPlAngV0CCor2 = TMath::ATan2(QycCor2, QxcCor2)/2.;
1183  evPlAngV0ACor3 = TMath::ATan2(QyaCor3, QxaCor3)/3.;
1184  evPlAngV0CCor3 = TMath::ATan2(QycCor3, QxcCor3)/3.;
1185  }
1186  else{
1187  evPlAngV0ACor2 = TMath::ATan2(Qya2, Qxa2)/2.;
1188  evPlAngV0CCor2 = TMath::ATan2(Qyc2, Qxc2)/2.;
1189  evPlAngV0ACor3 = TMath::ATan2(Qya3, Qxa3)/3.;
1190  evPlAngV0CCor3 = TMath::ATan2(Qyc3, Qxc3)/3.;
1191  }
1192  }
1193  else{
1194  AliEventplane *ep = aodEvent->GetEventplane();
1195  evPlAngV0ACor2 = ep->GetEventplane("V0A", aodEvent, 2);
1196  evPlAngV0CCor2 = ep->GetEventplane("V0C", aodEvent, 2);
1197  evPlAngV0ACor3 = ep->GetEventplane("V0A", aodEvent, 3);
1198  evPlAngV0CCor3 = ep->GetEventplane("V0C", aodEvent, 3);
1199  }
1200 
1201 
1206 
1207  fHKsPhiEP->Fill(fZvtx,fgPsi2v0c);
1208  //loop track and get pid
1209  for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
1210  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(iT));
1211  if(!aodTrack) AliFatal("Not a standard AOD");
1212 
1213  if (!aodTrack){
1214  continue;
1215  }
1216 
1217  Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
1218  if(fFillDCA)
1219  trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
1220 
1221  if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
1222  continue;
1223  }
1224 
1225  Double_t b[2] = {-99., -99.};
1226  Double_t bCov[3] = {-99., -99., -99.};
1227 
1228  AliAODTrack param(*aodTrack);
1229  if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
1230  continue;
1231  }
1232 
1233  if (!fFillDCA && ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4)))
1234  continue;
1235 
1236  if(fFillDCA && (TMath::Abs(b[0]) > 3.0 || TMath::Abs(b[1]) > 3))
1237  continue;
1238 
1239  // re-map the container in an array to do the analysis for V0A and V0C within a loop
1240  Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1242  AliFlowVZEROQA *QA[2] = {fQA,fQA2};
1243 
1244  Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1246  AliFlowVZEROQA *QAv3[2] = {fQAv3,fQA2v3};
1247 
1248  // Fill MC results
1249  if(fIsMC && mcArray){
1250  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1251  Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1252 
1253  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
1254 
1255  Float_t v2mc = TMath::Cos(2*(aodTrack->Phi() - evplaneMC));
1256 
1257  fContAllChargesMC->Fill(0,aodTrack->Pt(),v2mc,xMC);
1258 
1259  Int_t iS = TMath::Abs(((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->GetPdgCode());
1260  if(iS==11){
1261  fContAllChargesMC->Fill(4,aodTrack->Pt(),v2mc,xMC);
1262  }
1263  else if(iS==13){
1264  fContAllChargesMC->Fill(5,aodTrack->Pt(),v2mc,xMC);
1265  }
1266  else if(iS==211){
1267  fContAllChargesMC->Fill(1,aodTrack->Pt(),v2mc,xMC);
1268  }
1269  else if(iS==321){
1270  fContAllChargesMC->Fill(2,aodTrack->Pt(),v2mc,xMC);
1271  }
1272  else if(iS==2212){
1273  fContAllChargesMC->Fill(3,aodTrack->Pt(),v2mc,xMC);
1274  }
1275  }
1276 
1277  for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1278 
1279  if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
1280 
1281  Float_t v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1282  Float_t v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1283 
1284  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1285  Float_t dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1286  Float_t *probRead = fPID->GetProb();
1287  Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1288  Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1289  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
1290  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
1291 
1292  // in case fill DCA info
1293  if(fFillDCA){
1294  if(TMath::Abs(b[0]) > 0.1){
1295  x[5] = 1;
1296  x3[5] = 1;
1297  }
1298  if(TMath::Abs(b[0]) > 0.3){
1299  x[5] = 2;
1300  x3[5] = 2;
1301  }
1302  if(fIsMC && mcArray){
1303  if(!((AliAODMCParticle*)mcArray->At(TMath::Abs(aodTrack->GetLabel())))->IsPhysicalPrimary()){
1304  x[5] += 3;
1305  x3[5] += 3;
1306  }
1307  }
1308  }
1309 
1310  // Fill no PID
1311  if(fV2) contV0[iV0]->Fill(0,aodTrack->Pt(),v2V0,x);
1312  if(fV3) contV0v3[iV0]->Fill(0,aodTrack->Pt(),v3V0,x3);
1313 
1314 
1315  Double_t dedxExp[8];
1316  Float_t tof = -1;
1317  Double_t inttimes[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1318  Double_t expTOFsigma[8] = {0.,0.,0.,0.,0.,0.,0.,0.};
1319 
1320  Float_t nsigmaTPC[8];
1321  Float_t nsigmaTOF[8];
1322 
1323  if(aodTrack->GetDetPid()){ // check the PID object is available
1324  for(Int_t iS=0;iS < 8;iS++){
1325  dedxExp[iS] = fPID->GetExpDeDx(aodTrack,iS);
1326  nsigmaTPC[iS] = (dedx - dedxExp[iS])/(dedxExp[iS]*0.07);
1327  // printf("TPC %i = %f (%f %f)\n",iS, nsigmaTPC[iS],dedx, dedxExp[iS]);
1328  }
1329 
1330  if(fPID->GetCurrentMask(1)){ // if TOF is present
1331  Float_t ptrack = aodTrack->P();
1332  tof = aodTrack->GetTOFsignal() - fPID->GetESDpid()->GetTOFResponse().GetStartTime(ptrack);
1333  aodTrack->GetIntegratedTimes(inttimes);
1334 
1335  for(Int_t iS=5;iS < 8;iS++) // extra info for light nuclei
1336  inttimes[iS] = inttimes[0] / ptrack * mass[iS] * TMath::Sqrt(1+ptrack*ptrack/mass[iS]/mass[iS]);
1337 
1338  for(Int_t iS=0;iS<8;iS++){
1339  expTOFsigma[iS] = fPID->GetESDpid()->GetTOFResponse().GetExpectedSigma(ptrack, inttimes[iS], mass[iS]);
1340  nsigmaTOF[iS] = (tof - inttimes[iS])/expTOFsigma[iS];
1341  // printf("TOF %i = %f\n",iS, nsigmaTOF[iS]);
1342  }
1343  }
1344  }
1345 
1346  Float_t deltaPhiV0 = aodTrack->Phi() - evPlAngV0[iV0];
1347  if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1348  else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1349  if(deltaPhiV0 > TMath::Pi()) deltaPhiV0 -= 2*TMath::Pi();
1350  else if(deltaPhiV0 < -TMath::Pi()) deltaPhiV0 += 2*TMath::Pi();
1351 
1352  Float_t deltaPhiV0v3 = aodTrack->Phi() - evPlAngV0v3[iV0];
1353  if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1354  else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1355  if(deltaPhiV0v3 > TMath::Pi()) deltaPhiV0v3 -= 2*TMath::Pi();
1356  else if(deltaPhiV0v3 < -TMath::Pi()) deltaPhiV0v3 += 2*TMath::Pi();
1357 
1358  // variable to fill QA container
1359  Float_t xQA[5] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt()), 0.0,deltaPhiV0,x[4]}; // v2
1360  Float_t xQA3[5] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt()), 0.0,deltaPhiV0v3,x[4]}; // v3
1361 
1362  // extra QA TProfiles
1363  if(iV0==1 && aodTrack->Pt() < 20 && fPID->GetCurrentMask(0) && fPID->GetCurrentMask(1)){
1364  Float_t xQApid[2] = {static_cast<Float_t>(iC),static_cast<Float_t>(aodTrack->Pt())};
1365  fContQApid->Fill(0,nsigmaTPC[2],v2V0,xQApid); // v2 TPC (V0C) w.r.t pions
1366  fContQApid->Fill(1,nsigmaTOF[2],v2V0,xQApid); // v2 TOF (V0C) w.r.t. pions
1367  fContQApid->Fill(2,nsigmaTPC[3],v2V0,xQApid); // v2 TPC (V0C) w.r.t kaons
1368  fContQApid->Fill(3,nsigmaTOF[3],v2V0,xQApid); // v2 TOF (V0C) w.r.t. kaons
1369  fContQApid->Fill(4,nsigmaTPC[4],v2V0,xQApid); // v2 TPC (V0C) w.r.t protons
1370  fContQApid->Fill(5,nsigmaTOF[4],v2V0,xQApid); // v2 TOF (V0C) w.r.t. protons
1371  }
1372 
1373  // QA fill
1374  if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid() || dedx < 10. || aodTrack->Pt() < 0 || aodTrack->Pt() > 7){}
1375  else{
1376  if(TMath::Abs(nsigmaTPC[2])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[2])<5))){ //pi
1377  xQA[2] = prob[2];
1378  xQA3[2] = xQA[2];
1379  if(fQAsw && fV2) QA[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA);
1380  if(fQAsw && fV3) QAv3[iV0]->Fill(0,nsigmaTPC[2],nsigmaTOF[2],xQA3);
1381  }
1382  if(TMath::Abs(nsigmaTPC[3])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[3])<5))){ //K
1383  xQA[2] = prob[3];
1384  xQA3[2] = xQA[2];
1385  if(fQAsw && fV2) QA[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA);
1386 // if(fQAsw && fV3) QAv3[iV0]->Fill(1,nsigmaTPC[3],nsigmaTOF[3],xQA3);
1387  }
1388  if(TMath::Abs(nsigmaTPC[4])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[4])<5))){//p
1389  xQA[2] = prob[4];
1390  xQA3[2] = xQA[2];
1391  if(fQAsw && fV2) QA[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA);
1392 // if(fQAsw && fV3) QAv3[iV0]->Fill(2,nsigmaTPC[4],nsigmaTOF[4],xQA3);
1393  }
1394  if(TMath::Abs(nsigmaTPC[0])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[0])<5))){//e
1395  xQA[2] = prob[0];
1396  xQA3[2] = xQA[2];
1397 // if(fQAsw && fV2) QA[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA);
1398 // if(fQAsw && fV3) QAv3[iV0]->Fill(3,nsigmaTPC[0],nsigmaTOF[0],xQA3);
1399  }
1400  if(TMath::Abs(nsigmaTPC[5])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[5])<5))){//d
1401  xQA[2] = prob[5];
1402  xQA3[2] = xQA[2];
1403  // if(fQAsw && fV2) QA[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA);
1404  // if(fQAsw && fV3) QAv3[iV0]->Fill(4,nsigmaTPC[5],nsigmaTOF[5],xQA3);
1405  }
1406  if(TMath::Abs(nsigmaTPC[6])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[6])<5))){//t
1407  xQA[2] = prob[6];
1408  xQA3[2] = xQA[2];
1409  // if(fQAsw && fV2) QA[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA);
1410  // if(fQAsw && fV3) QAv3[iV0]->Fill(5,nsigmaTPC[6],nsigmaTOF[6],xQA3);
1411  }
1412  if(TMath::Abs(nsigmaTPC[7])<5 && (!(fPID->GetCurrentMask(1)) || (TMath::Abs(nsigmaTOF[7])<5))){//He3
1413  xQA[2] = prob[7];
1414  xQA3[2] = xQA[2];
1415  // if(fQAsw && fV2) QA[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA);
1416  // if(fQAsw && fV3) QAv3[iV0]->Fill(6,nsigmaTPC[7],nsigmaTOF[7],xQA3);
1417  }
1418  }
1419 
1420  //pid selection
1421  if(!(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID and PID object strictly required (very important!!!!)
1422  else if(prob[2] > 0.6){ // pi
1423  x[2] = prob[2];
1424  x3[2] = x[2];
1425  if(TMath::Abs(nsigmaTPC[2]) < 5){ // TPC 5 sigma extra cut to accept the track
1426  if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1427  if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1428  if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][0]->Fill(aodTrack->Pt(),b[0]);
1429  else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][0]->Fill(aodTrack->Pt(),b[0]);
1430  }
1431  }
1432  else if(prob[3] > 0.6){ // K
1433  x[2] = prob[3];
1434  x3[2] = x[2];
1435  if(TMath::Abs(nsigmaTPC[3]) < 5){
1436  if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1437  if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1438  if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][1]->Fill(aodTrack->Pt(),b[0]);
1439  else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][1]->Fill(aodTrack->Pt(),b[0]);
1440  }
1441  }
1442  else if(prob[4] > 0.6){ // p
1443  x[2] = prob[4];
1444  x3[2] = x[2];
1445  if(TMath::Abs(nsigmaTPC[4]) < 5){
1446  if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1447  if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1448  if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][2]->Fill(aodTrack->Pt(),b[0]);
1449  else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][2]->Fill(aodTrack->Pt(),b[0]);
1450  }
1451  }
1452  else if(prob[0] > 0.6){ // e
1453  x[2] = prob[0];
1454  x3[2] = x[2];
1455  if(TMath::Abs(nsigmaTPC[0]) < 5){
1456  if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1457  if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1458  if(x[2] > 0.9 && x[5] < 3) fHdcaPt[iC][3]->Fill(aodTrack->Pt(),b[0]);
1459  else if(x[2] > 0.9 && fIsMC) fHdcaPtSec[iC][3]->Fill(aodTrack->Pt(),b[0]);
1460  }
1461  }
1462  else if(prob[1] > 0.6){ // mu
1463  x[2] = prob[1];
1464  x3[2] = x[2];
1465  if(TMath::Abs(nsigmaTPC[1]) < 5){
1466  if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1467  if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1468  }
1469  }
1470  else if(prob[5] > 0.6){ // d
1471  x[2] = prob[5];
1472  x3[2] = x[2];
1473  if(TMath::Abs(nsigmaTPC[5]) < 5){
1474  if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1475  if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1476  if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][4]->Fill(aodTrack->Pt(),b[0]);
1477  else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][4]->Fill(aodTrack->Pt(),b[0]);
1478  }
1479  }
1480  else if(prob[6] > 0.6){ // t
1481  x[2] = prob[6];
1482  x3[2] = x[2];
1483  if(TMath::Abs(nsigmaTPC[6]) < 5){
1484  if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1485  if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1486  if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][5]->Fill(aodTrack->Pt(),b[0]);
1487  else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][5]->Fill(aodTrack->Pt(),b[0]);
1488  }
1489  }
1490  else if(prob[7] > 0.6){ // He3
1491  x[2] = prob[7];
1492  x3[2] = x[2];
1493  if(TMath::Abs(nsigmaTPC[7]) < 5){
1494  if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1495  if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1496  if(x[2] > 0.9 && x[5] < 3 && x[1] < 0) fHdcaPt[iC][6]->Fill(aodTrack->Pt(),b[0]);
1497  else if(x[2] > 0.9 && fIsMC && x[1] < 0) fHdcaPtSec[iC][6]->Fill(aodTrack->Pt(),b[0]);
1498  }
1499  }
1500 
1501  if(x[4] > 0.5){ // if TOF was present redo TPC stand alone PID to check the PID in the same acceptance (PID mask = 2)
1502  fPID->ResetDetOR(1); // exclude TOF from PID
1503  tofMismProb = 0;
1504 
1505  fPID->ComputeProb(aodTrack,fOutputAOD);
1506  dedx = fPID->GetDeDx();//aodTrack->GetTPCsignal();
1507  probRead = fPID->GetProb();
1508 
1509  fPID->SetDetOR(1); // include TOF for PID
1510  }
1511  Float_t probTPC[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]}; // TPC stand alone prbabilities
1512 
1513  //pid selection TPC S.A. with TOF matching
1514  x[4]*=2; // set the mask to 2 id TOF is present
1515  if(x[4]<1 || !(fPID->GetCurrentMask(0)) || !aodTrack->GetDetPid()){} // TPC PID S.A. PID in TOF acceptance
1516  else if(probTPC[2] > 0.6){ // pi
1517  x[2] = probTPC[2];
1518  x3[2] = x[2];
1519  if(TMath::Abs(nsigmaTPC[2]) < 5){
1520  if(fV2) contV0[iV0]->Fill(1,aodTrack->Pt(),v2V0,x);
1521  if(fV3) contV0v3[iV0]->Fill(1,aodTrack->Pt(),v3V0,x3);
1522  }
1523  }
1524  else if(probTPC[3] > 0.6){ // K
1525  x[2] = probTPC[3];
1526  x3[2] = x[2];
1527  if(TMath::Abs(nsigmaTPC[3]) < 5){
1528  if(fV2) contV0[iV0]->Fill(2,aodTrack->Pt(),v2V0,x);
1529  if(fV3) contV0v3[iV0]->Fill(2,aodTrack->Pt(),v3V0,x3);
1530  }
1531  }
1532  else if(probTPC[4] > 0.6){ // p
1533  x[2] = probTPC[4];
1534  x3[2] = x[2];
1535  if(TMath::Abs(nsigmaTPC[4]) < 5){
1536  if(fV2) contV0[iV0]->Fill(3,aodTrack->Pt(),v2V0,x);
1537  if(fV3) contV0v3[iV0]->Fill(3,aodTrack->Pt(),v3V0,x3);
1538  }
1539  }
1540  else if(probTPC[0] > 0.6){ // e
1541  x[2] = probTPC[0];
1542  x3[2] = x[2];
1543  if(TMath::Abs(nsigmaTPC[0]) < 5){
1544  if(fV2) contV0[iV0]->Fill(4,aodTrack->Pt(),v2V0,x);
1545  if(fV3) contV0v3[iV0]->Fill(4,aodTrack->Pt(),v3V0,x3);
1546  }
1547  }
1548  else if(probTPC[1] > 0.6){ // mu
1549  x[2] = probTPC[1];
1550  x3[2] = x[2];
1551  if(TMath::Abs(nsigmaTPC[1]) < 5){
1552  if(fV2) contV0[iV0]->Fill(8,aodTrack->Pt(),v2V0,x);
1553  if(fV3) contV0v3[iV0]->Fill(8,aodTrack->Pt(),v3V0,x3);
1554  }
1555  }
1556  else if(probTPC[5] > 0.6){ // d
1557  x[2] = probTPC[5];
1558  x3[2] = x[2];
1559  if(TMath::Abs(nsigmaTPC[5]) < 5){
1560  if(fV2) contV0[iV0]->Fill(5,aodTrack->Pt(),v2V0,x);
1561  if(fV3) contV0v3[iV0]->Fill(5,aodTrack->Pt(),v3V0,x3);
1562  }
1563  }
1564  else if(probTPC[6] > 0.6){ // t
1565  x[2] = probTPC[6];
1566  x3[2] = x[2];
1567  if(TMath::Abs(nsigmaTPC[6]) < 5){
1568  if(fV2) contV0[iV0]->Fill(6,aodTrack->Pt(),v2V0,x);
1569  if(fV3) contV0v3[iV0]->Fill(6,aodTrack->Pt(),v3V0,x3);
1570  }
1571  }
1572  else if(probTPC[7] > 0.6){ // He3
1573  x[2] = probTPC[7];
1574  x3[2] = x[2];
1575  if(TMath::Abs(nsigmaTPC[7]) < 5){
1576  if(fV2) contV0[iV0]->Fill(7,aodTrack->Pt()*2,v2V0,x);
1577  if(fV3) contV0v3[iV0]->Fill(7,aodTrack->Pt()*2,v3V0,x3);
1578  }
1579  }
1580  } // end side loop
1581  } // end track loop
1582 
1583  // my V0 loop
1584  for(Int_t imy=0;imy<fNK0s;imy++){
1585  Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1586  Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1587 
1590 
1591  for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1592  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
1593  Float_t x3[6] = {static_cast<Float_t>(iC),-1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
1594 
1595  Float_t v2V0 = TMath::Cos(2*(fPhiK0s[imy] - evPlAngV0[iV0]));
1596  Float_t v3V0 = TMath::Cos(3*(fPhiK0s[imy] - evPlAngV0v3[iV0]));
1597  if(fV2) contV0[iV0]->Fill(9,fPtK0s[imy],v2V0,x);
1598  if(fV3) contV0v3[iV0]->Fill(9,fPtK0s[imy],v3V0,x3);
1599  }
1600  }
1601 
1602  // V0 loop
1603  Int_t nV0s = fOutputAOD->GetNumberOfV0s();
1604  AliAODv0 *myV0;
1605  Double_t /*dQT, dALPHA, dPT,*/ dMASS=0.0;
1606  for (Int_t i=0; i!=nV0s; ++i) {
1607  myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
1608  if(!myV0) continue;
1609  if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > fEtaCut) continue; // skipping low momentum
1610  Int_t pass = PassesAODCuts(myV0,fOutputAOD,0);
1611  if(pass) {
1612  dMASS = myV0->MassK0Short();
1613  pass = 3;
1614  fHK0sMass2->Fill(myV0->Pt(),dMASS);
1615  }
1616  if(TMath::Abs(dMASS-0.497)/0.005 > 3){
1617  pass = PassesAODCuts(myV0,fOutputAOD,1);
1618  if(pass) dMASS = myV0->MassLambda();
1619  if(pass==2) dMASS = myV0->MassAntiLambda();
1620  }
1621 
1622  if(pass){// 1 lambda, 2 antilambda, 3=K0s
1623  //dPT=myV0->Pt();
1624  //dQT=myV0->PtArmV0();
1625  //dALPHA=myV0->AlphaV0();
1626 
1627  Int_t iPos, iNeg;
1628  AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0);
1629  if(iT->Charge()>0) {
1630  iPos = 0; iNeg = 1;
1631  } else {
1632  iPos = 1; iNeg = 0;
1633  }
1634 
1635  // check if one of the daugthers was already used
1636  if(pass == 3 && TMath::Abs(dMASS-0.497)/0.005 < 1){
1637  fHKsPhi->Fill(fZvtx, myV0->Phi());
1638  }
1639 
1640  if(pass == 1000){ // disable
1641  Bool_t used = kFALSE;
1642  for(Int_t ii=0;ii<nusedForK0s;ii++){
1643  if(myV0->GetDaughter(iNeg) == usedForK0s1[ii] || myV0->GetDaughter(iPos) == usedForK0s2[ii]){
1644  used = kTRUE;
1645  }
1646  }
1647  if((!used) && nusedForK0s < 1000){
1648  nusedForK0s++;
1649  usedForK0s1[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iNeg);
1650  usedForK0s2[nusedForK0s] = (AliAODTrack *) myV0->GetDaughter(iPos);
1651  printf("accepted\n");
1652  }
1653  else{
1654  dMASS = 0;
1655  printf("rejected\n");
1656  }
1657  }
1658 
1659  iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1660  AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1661 
1662  // re-map the container in an array to do the analysis for V0A and V0C within a loop
1663  Float_t evPlAngV0[2] = {evPlAngV0ACor2,evPlAngV0CCor2};
1665 
1666  Float_t evPlAngV0v3[2] = {evPlAngV0ACor3,evPlAngV0CCor3};
1668 
1669  for(Int_t iV0=0;iV0<2;iV0++){ // loop on A and C side
1670 
1671  if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAngV0[iV0],evPlRes[iV0*8+iC]); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
1672 
1673  Float_t v2V0 = TMath::Cos(2*(myV0->Phi() - evPlAngV0[iV0]));
1674  Float_t v3V0 = TMath::Cos(3*(myV0->Phi() - evPlAngV0v3[iV0]));
1675 
1676  Float_t deltaphi = myV0->Phi()- evPlAngV0[iV0];
1677  if(deltaphi > TMath::Pi()) deltaphi -= 2*TMath::Pi();
1678  if(deltaphi < -TMath::Pi()) deltaphi += 2*TMath::Pi();
1679 
1680  Float_t x[6] = {static_cast<Float_t>(iC),1,1,evPlAngV0[iV0],1,0}; // to fill analysis v2 container
1681  Float_t x3[6] = {static_cast<Float_t>(iC),1,1,evPlAngV0v3[iV0],1,0}; // to fill analysis v3 container
1682 
1683  Float_t decaylength = myV0->DecayLengthXY(fOutputAOD->GetPrimaryVertex());
1684  // printf("decay length = %f\n",decaylength);
1685 
1686  if(pass==2){ // anti-lambda charge = -1
1687  x[1] = -1;
1688  x3[1] = -1;
1689  }
1690 
1691  if(decaylength < fMinDistV0) pass = 0;
1692  if(decaylength > fMaxDistV0) pass = 0;
1693 
1694  Float_t nsigma = 0;
1695  if(pass < 3)
1696  nsigma = TMath::Abs(dMASS-1.116)/0.0016;
1697  else if(pass == 3)
1698  nsigma = TMath::Abs(dMASS-0.497)/0.005;
1699 
1700  if(nsigma < 1)
1701  x[2] = 0.95;
1702  else if(nsigma < 2)
1703  x[2] = 0.85;
1704  else if(nsigma < 3)
1705  x[2] = 0.75;
1706  else if(nsigma < 4)
1707  x[2] = 0.65;
1708  else
1709  x[2] = 0.5;
1710 
1711  x3[2] = x[2];
1712 
1713  // Fill Container for lambda and Ks
1714  if(fV2 && pass == 3 && x[2] > 0.6){
1715  contV0[iV0]->Fill(9,myV0->Pt(),v2V0,x);
1716  fHctauPtEP->Fill(myV0->Pt(),deltaphi,decaylength);//ciao
1717  if(myV0->Pt() < 1.1 && myV0->Pt() > 0.9) fHctauAt1EP->Fill(decaylength,deltaphi);
1718  }
1719  if(fV3 && pass == 3 && x[2] > 0.6) contV0v3[iV0]->Fill(9,myV0->Pt(),v3V0,x3);
1720  if(fV2 && pass < 3 && x[2] > 0.6) contV0[iV0]->Fill(10,myV0->Pt(),v2V0,x);
1721  if(fV3 && pass < 3 && x[2] > 0.6) contV0v3[iV0]->Fill(10,myV0->Pt(),v3V0,x3);
1722 
1723  if(pass < 3){ // lambda
1724  AliAODTrack* aodTrack = iT;
1725  if(pass==2) aodTrack=jT;
1726 
1727  v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1728  v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1729 
1730  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1731  Float_t *probRead = fPID->GetProb();
1732  Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1733  Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1734 
1735  if(prob[4] < 0.61) prob[4] = 0.61;
1736 
1737  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
1738  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
1739 
1740  // Fill Container for (anti)proton from lambda
1741  if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1742  if(fV2) contV0[iV0]->Fill(11,aodTrack->Pt(),v2V0,xdec);
1743  if(fV3) contV0v3[iV0]->Fill(11,aodTrack->Pt(),v3V0,xdec3);
1744  }
1745  }
1746  else if(pass == 3){
1747  AliAODTrack* aodTrack = iT;
1748 
1749  v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1750  v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1751 
1752  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1753  Float_t *probRead = fPID->GetProb();
1754  Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
1755  Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1756 
1757  if(prob[2] < 0.61) prob[2] = 0.61;
1758 
1759  Float_t xdec[6] = {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
1760  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
1761 
1762  if(nsigma < 2 && xdec[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1763  if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdec);
1764  if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdec3);
1765  }
1766 
1767  aodTrack = jT;
1768  v2V0 = TMath::Cos(2*(aodTrack->Phi() - evPlAngV0[iV0]));
1769  v3V0 = TMath::Cos(3*(aodTrack->Phi() - evPlAngV0v3[iV0]));
1770 
1771  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
1772  Float_t *probRead2 = fPID->GetProb();
1773  Float_t prob2[8] = {probRead2[0],probRead2[1],probRead2[2],probRead2[3],probRead2[4],probRead2[5],probRead2[6],probRead2[7]};
1774  Float_t tofMismProb2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
1775 
1776  if(prob2[2] < 0.61) prob2[2] = 0.61;
1777 
1778  Float_t xdecB[6] = {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
1779  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
1780 
1781  if(nsigma < 2 && xdecB[2] > 0.6 && TMath::Abs(aodTrack->Eta()) < 0.8){
1782  if(fV2) contV0[iV0]->Fill(12,aodTrack->Pt(),v2V0,xdecB);
1783  if(fV3) contV0v3[iV0]->Fill(12,aodTrack->Pt(),v3V0,xdecB3);
1784  }
1785  }
1786  }
1787 
1788  }
1789  } // end loop on V0
1790 
1791 
1792  // Fill EP distribution histograms
1793  if(fV2) fPhiRPv0A->Fill(iC,evPlAngV0ACor2);
1794  if(fV2) fPhiRPv0C->Fill(iC,evPlAngV0CCor2);
1795 
1796  if(fV3) fPhiRPv0Av3->Fill(iC,evPlAngV0ACor3);
1797  if(fV3) fPhiRPv0Cv3->Fill(iC,evPlAngV0CCor3);
1798 
1799  // Fill histograms needed for resolution evaluation
1800  if(fV2) fHResTPCv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0ACor2)));
1801  if(fV2) fHResTPCv0C2->Fill(Double_t(iC), TMath::Cos(2*(evPlAng2 - evPlAngV0CCor2)));
1802  if(fV2) fHResv0Cv0A2->Fill(Double_t(iC), TMath::Cos(2*(evPlAngV0ACor2 - evPlAngV0CCor2)));
1803 
1804  if(fV3) fHResTPCv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0ACor3)));
1805  if(fV3) fHResTPCv0C3->Fill(Double_t(iC), TMath::Cos(3*(evPlAng3 - evPlAngV0CCor3)));
1806  if(fV3) fHResv0Cv0A3->Fill(Double_t(iC), TMath::Cos(3*(evPlAngV0ACor3 - evPlAngV0CCor3)));
1807  }
1808 
1809 
1810 
1811  // clean track array
1812  for(Int_t i=0;i < nusedForK0s;i++){
1813  usedForK0s1[i] = NULL;
1814  usedForK0s2[i] = NULL;
1815  }
1816 }
1817 
1818 //_____________________________________________________________________________
1820 {
1821 
1822  Float_t zvtx = -999;
1823 
1824  const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
1825  if (!vtxAOD)
1826  return zvtx;
1827  if(vtxAOD->GetNContributors()>0)
1828  zvtx = vtxAOD->GetZ();
1829 
1830  return zvtx;
1831 }
1832 //_____________________________________________________________________________
1834 {
1835  // Terminate loop
1836  Printf("Terminate()");
1837 }
1838 //_____________________________________________________________________________
1840  TString oadbfilename = "$ALICE_PHYSICS/OADB/PWGCF/VZERO/VZEROcalibEP.root";
1841  TFile *foadb = TFile::Open(oadbfilename.Data());
1842 
1843  if(!foadb){
1844  printf("OADB file %s cannot be opened\n",oadbfilename.Data());
1845  return;
1846  }
1847 
1848  AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
1849  if(!cont){
1850  printf("OADB object hMultV0BefCorr is not available in the file\n");
1851  return;
1852  }
1853 
1854  if(!(cont->GetObject(run))){
1855  printf("OADB object hMultV0BefCorr is not available for run %i (used run 137366)\n",run);
1856  run = 137366;
1857  }
1858  fMultV0 = ((TH2F *) cont->GetObject(run))->ProfileX();
1859 
1860  TF1 *fpol0 = new TF1("fpol0","pol0");
1861  fMultV0->Fit(fpol0,"","",0,31);
1862  fV0Cpol = fpol0->GetParameter(0);
1863  fMultV0->Fit(fpol0,"","",32,64);
1864  fV0Apol = fpol0->GetParameter(0);
1865 
1866  for(Int_t iside=0;iside<2;iside++){
1867  for(Int_t icoord=0;icoord<2;icoord++){
1868  for(Int_t i=0;i < 9;i++){
1869  char namecont[100];
1870  if(iside==0 && icoord==0)
1871  snprintf(namecont,100,"hQxc2_%i",i);
1872  else if(iside==1 && icoord==0)
1873  snprintf(namecont,100,"hQxa2_%i",i);
1874  else if(iside==0 && icoord==1)
1875  snprintf(namecont,100,"hQyc2_%i",i);
1876  else if(iside==1 && icoord==1)
1877  snprintf(namecont,100,"hQya2_%i",i);
1878 
1879  cont = (AliOADBContainer*) foadb->Get(namecont);
1880  if(!cont){
1881  printf("OADB object %s is not available in the file\n",namecont);
1882  return;
1883  }
1884 
1885  if(!(cont->GetObject(run))){
1886  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1887  run = 137366;
1888  }
1889  fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1890  fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1891 
1892  //for v3
1893  if(iside==0 && icoord==0)
1894  snprintf(namecont,100,"hQxc3_%i",i);
1895  else if(iside==1 && icoord==0)
1896  snprintf(namecont,100,"hQxa3_%i",i);
1897  else if(iside==0 && icoord==1)
1898  snprintf(namecont,100,"hQyc3_%i",i);
1899  else if(iside==1 && icoord==1)
1900  snprintf(namecont,100,"hQya3_%i",i);
1901 
1902  cont = (AliOADBContainer*) foadb->Get(namecont);
1903  if(!cont){
1904  printf("OADB object %s is not available in the file\n",namecont);
1905  return;
1906  }
1907 
1908  if(!(cont->GetObject(run))){
1909  printf("OADB object %s is not available for run %i (used run 137366)\n",namecont,run);
1910  run = 137366;
1911  }
1912  fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
1913  fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
1914 
1915  }
1916  }
1917  }
1918 }
1919 //=======================================================================
1921 {
1922  Int_t set = 0;
1923  Float_t fV0Cuts[9];
1924  // defines cuts to be used
1925  // fV0Cuts[9] dl dca ctp d0 d0d0 qt minEta maxEta PID
1926  switch(set) {
1927  case(0): // No cuts
1928  fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1929  fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1930  fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 0;
1931  break;
1932  case(1): // Tight cuts
1933  fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1934  fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1935  fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 0;
1936  break;
1937  case(2): // Tight cuts + PID
1938  fV0Cuts[0] = +0.5; fV0Cuts[1] = +0.5; fV0Cuts[2] = +0.998;
1939  fV0Cuts[3] = +0.1; fV0Cuts[4] = +0.0; fV0Cuts[5] = +0.105;
1940  fV0Cuts[6] = -0.8; fV0Cuts[7] = +0.8; fV0Cuts[8] = 1;
1941  break;
1942  case(3): // No cuts + PID
1943  fV0Cuts[0] = -1e+6; fV0Cuts[1] = +1e+6; fV0Cuts[2] = -1e+6;
1944  fV0Cuts[3] = -1e+6; fV0Cuts[4] = +1e+6; fV0Cuts[5] = -1e+6;
1945  fV0Cuts[6] = -1e+6; fV0Cuts[7] = +1e+6; fV0Cuts[8] = 1;
1946  break;
1947  }
1948 
1949  // daughter cuts
1950  if(! fCutsDaughter){
1951  fCutsDaughter = new AliESDtrackCuts(Form("daughter_cuts_%s","ESD") );
1952  fCutsDaughter->SetPtRange(0.2,10.0);
1953  fCutsDaughter->SetEtaRange(-0.8, 0.8 );
1954  fCutsDaughter->SetMinNClustersTPC(80);
1955  fCutsDaughter->SetMaxChi2PerClusterTPC(4.0);
1956  fCutsDaughter->SetRequireTPCRefit(kTRUE);
1957  fCutsDaughter->SetAcceptKinkDaughters(kFALSE);
1958  }
1959 
1960  if (myV0->GetOnFlyStatus() ) return 0;
1961  //the following is needed in order to evualuate track-quality
1962  AliAODTrack *iT, *jT;
1963  AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1964  Double_t pos[3],cov[6];
1965  vV0s->GetXYZ(pos);
1966  vV0s->GetCovarianceMatrix(cov);
1967  const AliESDVertex vESD(pos,cov,100.,100);
1968  // TESTING CHARGE
1969  int iPos, iNeg;
1970  iT=(AliAODTrack*) myV0->GetDaughter(0);
1971  if(iT->Charge()>0) {
1972  iPos = 0; iNeg = 1;
1973  } else {
1974  iPos = 1; iNeg = 0;
1975  }
1976  // END OF TEST
1977 
1978  iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1979  AliESDtrack ieT( iT );
1980  ieT.SetTPCClusterMap( iT->GetTPCClusterMap() );
1981  ieT.SetTPCSharedMap( iT->GetTPCSharedMap() );
1982  ieT.SetTPCPointsF( iT->GetTPCNclsF() );
1983  ieT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1984  if (!fCutsDaughter->IsSelected( &ieT ) ) return 0;
1985 
1986  jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1987  AliESDtrack jeT( jT );
1988  jeT.SetTPCClusterMap( jT->GetTPCClusterMap() );
1989  jeT.SetTPCSharedMap( jT->GetTPCSharedMap() );
1990  jeT.SetTPCPointsF( jT->GetTPCNclsF() );
1991  jeT.RelateToVertex(&vESD, tAOD->GetMagneticField(), 100);
1992  if (!fCutsDaughter->IsSelected( &jeT ) ) return 0;
1993 
1994  Double_t pvertex[3];
1995  pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1996  pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1997  pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1998  Double_t dDL=myV0->DecayLengthV0( pvertex );
1999  Double_t dDCA=myV0->DcaV0Daughters();
2000  Double_t dCTP=myV0->CosPointingAngle( pvertex );
2001  Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2002  Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
2003  Double_t dD0D0=dD0P*dD0M;
2004  Double_t dQT=myV0->PtArmV0();
2005  Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
2006  if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
2007 // Double_t dPT=myV0->Pt();
2008  Double_t dETA=myV0->Eta();
2009  Int_t passes = 1;
2010  if(dDL <fV0Cuts[0]) passes = 0;
2011  if(dDCA >fV0Cuts[1]) passes = 0;
2012  if(dCTP <fV0Cuts[2]) passes = 0;
2013  if(TMath::Abs(dD0P) <fV0Cuts[3]) passes = 0;
2014  if(TMath::Abs(dD0M) <fV0Cuts[3]) passes = 0;
2015  if(dD0D0>fV0Cuts[4]) passes = 0;
2016  if(dETA <fV0Cuts[6]) passes = 0;
2017  if(dETA >fV0Cuts[7]) passes = 0;
2018  if(specie==0) if(dQT<fV0Cuts[5]) passes = 0;
2019  if(specie==1&&passes==1&&dALPHA<0) passes = 2; // antilambda
2020 
2021 
2022 // if(jT->Pt() < 0.5*myV0->Pt() || iT->Pt() < 0.5*myV0->Pt()) passes = 0;
2023 
2024 
2025  // additional cut
2026 // if(!(iT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2027 // if(!(jT->GetStatus() & AliAODTrack::kTPCrefit)) passes = 0;
2028 
2029 // if(!(iT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2030 // if(!(jT->GetStatus() & AliAODTrack::kITSrefit)) passes = 0;
2031 
2032 // if(!(iT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2033 // if(!(jT->GetStatus() & AliAODTrack::kTOFout)) passes = 0;
2034 
2035  Bool_t trkFlag = iT->TestFilterBit(1); // TPC only tracks (4,global track)
2036  Bool_t trkFlag2 = jT->TestFilterBit(1); // TPC only tracks (4,global track)
2037 
2038  if(!trkFlag) passes = 0;
2039  if(!trkFlag2) passes = 0;
2040 
2041  if(passes&&fV0Cuts[8]) {
2042 
2043  Double_t dedxExp[8];
2044  fPID->ComputeProb(iT,tAOD); // compute Bayesian probabilities
2045  Float_t nsigmaTPC[8];
2046 
2047 // Int_t tofMatch=0;
2048 // Int_t tofMatch2=0;
2049 
2050  if(iT->GetDetPid()){ // check the PID object is available
2051  for(Int_t iS=0;iS < 8;iS++){
2052  dedxExp[iS] = fPID->GetExpDeDx(iT,iS);
2053  nsigmaTPC[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2054  }
2055  }
2056  else{
2057  for(Int_t iS=0;iS < 8;iS++)
2058  nsigmaTPC[iS] = 10;
2059  }
2060 
2061 
2062 // if(fPID->GetCurrentMask(1)) // if TOF is present
2063 // tofMatch = 1;
2064 
2065 // Float_t tofMismProbMC = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2066 
2067  Float_t *probRead = fPID->GetProb();
2068  Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2069 
2070  fPID->ComputeProb(jT,tAOD); // compute Bayesian probabilities
2071  Float_t nsigmaTPC2[8];
2072  if(jT->GetDetPid()){ // check the PID object is available
2073  for(Int_t iS=0;iS < 8;iS++){
2074  dedxExp[iS] = fPID->GetExpDeDx(jT,iS);
2075  nsigmaTPC2[iS] = (fPID->GetDeDx() - dedxExp[iS])/(dedxExp[iS]*0.07);
2076  }
2077  }
2078  else{
2079  for(Int_t iS=0;iS < 8;iS++)
2080  nsigmaTPC2[iS] = 10;
2081  }
2082 
2083 // if(fPID->GetCurrentMask(1)) // if TOF is present
2084 // tofMatch2 = 1;
2085 
2086 // Float_t tofMismProbMC2 = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2087 
2088  probRead = fPID->GetProb();
2089  Float_t prob2[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2090 
2091  if(jT->GetTPCNcls() < fNcluster) passes = 0;
2092  else if(iT->GetTPCNcls() < fNcluster) passes = 0;
2093 
2094 // if(! (tofMatch && tofMatch2)) passes = 0;
2095 
2096  /*
2097  Float_t dMASS = myV0->MassK0Short();
2098  Float_t nsigmaMass = TMath::Abs(dMASS-0.497)/0.005;
2099  if(specie == 0 && TMath::Abs(nsigmaMass) < 1 && myV0->Pt() > 1) printf("candidate i=(pt=%f-phi=%f-tof=%i) j=(pt=%f-phi=%f-tof=%i) \n",iT->Pt(),iT->Phi(),tofMatch,jT->Pt(),jT->Phi(),tofMatch2);
2100  */
2101 
2102  switch(specie) {
2103  case 0: // K0 PID
2104  if(0){
2105  if( ((jT->GetTPCmomentum()<15) &&
2106  (TMath::Abs(nsigmaTPC2[2])>3.)) || prob2[2] < 0.9)
2107  passes = 0;
2108  if( ((iT->GetTPCmomentum()<15) &&
2109  (TMath::Abs(nsigmaTPC[2])>3.))|| prob[2] < 0.9 )
2110  passes = 0;
2111  }
2112  break;
2113  case 1: // Lambda PID i==pos j ==neg
2114  if(passes==1) {
2115  if( (iT->GetTPCmomentum()<15) &&
2116  (TMath::Abs(nsigmaTPC[4])>3.) )
2117  passes = 0;
2118  if( (jT->GetTPCmomentum()<15) &&
2119  (TMath::Abs(nsigmaTPC2[2])>3.) )
2120  passes = 0;
2121  }
2122  if(passes==2) {
2123  if( (iT->GetTPCmomentum()<15) &&
2124  (TMath::Abs(nsigmaTPC[2])>3.) )
2125  passes = 0;
2126  if( (jT->GetTPCmomentum()<15) &&
2127  (TMath::Abs(nsigmaTPC2[4])>3.) )
2128  passes = 0;
2129  }
2130  break;
2131  }
2132  }
2133  return passes;
2134 }
2135 //=======================================================================
2137  fNK0s=0;
2138  fNpiPos=0;
2139  fNpiNeg=0;
2140 
2141  if(fModulationDEDx) fPID->SetPsiCorrectionDeDx(evPlAng2,1.0); // set the PID dE/dx correction as a function of the v2-EP (resolution is needed)
2142 
2143  // fill pion stacks
2144  Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
2145  for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
2146  AliAODTrack* aodTrack = dynamic_cast<AliAODTrack*>(fOutputAOD->GetTrack(iT));
2147  if(!aodTrack) AliFatal("Not a standard AOD");
2148 
2149  if (!aodTrack){
2150  continue;
2151  }
2152 
2153  Bool_t trkFlag = aodTrack->TestFilterBit(1); // TPC only tracks
2154 // trkFlag = aodTrack->TestFilterBit(4); // Global track, DCA loose cut
2155 
2156  if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
2157  continue;
2158  }
2159 
2160  Double_t b[2] = {-99., -99.};
2161  Double_t bCov[3] = {-99., -99., -99.};
2162 
2163  AliAODTrack param(*aodTrack);
2164  if (!param.PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov)){
2165  continue;
2166  }
2167 
2168  if(TMath::Abs(b[0]) < 0.5/aodTrack->Pt()) continue;
2169 
2170  fPID->ComputeProb(aodTrack,fOutputAOD); // compute Bayesian probabilities
2171  Float_t *probRead = fPID->GetProb();
2172  Float_t prob[8] = {probRead[0],probRead[1],probRead[2],probRead[3],probRead[4],probRead[5],probRead[6],probRead[7]};
2173  // Float_t tofMismProb = fPID->GetTOFMismProb(); // TOF mismatch probability requested to be lower than 50% for TOF analysis
2174 
2175  Int_t charge = aodTrack->Charge();
2176  if(prob[2] > 0.9){
2177  if(charge > 0){
2178  fIPiPos[fNpiPos] = iT;
2179  fNpiPos++;
2180  }
2181  else{
2182  fIPiNeg[fNpiNeg] = iT;
2183  fNpiNeg++;
2184  }
2185  }
2186  }
2187 
2188  for(Int_t i=0;i < fNpiPos;i++){
2189  AliAODTrack *pip = dynamic_cast<AliAODTrack*>(fOutputAOD->GetTrack(fIPiPos[i]));
2190  if(!pip) AliFatal("Not a standard AOD");
2191  AliESDtrack pipE(pip);
2192 
2193  for(Int_t j=0;j < fNpiNeg;j++){
2194  AliAODTrack *pin = dynamic_cast<AliAODTrack*>(fOutputAOD->GetTrack(fIPiNeg[j]));
2195  if(!pin) AliFatal("Not a standard AOD");
2196  AliESDtrack pinE(pin);
2197 
2198  Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
2199 
2200  Double_t pPos[3];
2201  Double_t pNeg[3];
2202  pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
2203  pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
2204 
2205  Float_t length = (xp+xn)*0.5;
2206 
2207  Float_t pxs = pPos[0] + pNeg[0];
2208  Float_t pys = pPos[1] + pNeg[1];
2209  Float_t pzs = pPos[2] + pNeg[2];
2210  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);
2211 
2212  Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
2213  Float_t phi = TMath::ATan2(pys,pxs);
2214  Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
2215 
2216  // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
2217 
2218  if(mindist < 0.2&& length > 1 && length < 25){
2219  fHK0sMass->Fill(pt,mass);
2220 
2221  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);
2222  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);
2223 
2224  Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
2225  Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
2226 
2227  fHK0vsLambda->Fill(mass,TMath::Min(massaL,massaAL));
2228 
2229  if(TMath::Abs(mass-0.497)/0.005 < 1 && massaL > 1.15 && massaAL > 1.15){
2230  fPhiK0s[fNK0s] = phi;
2231  fPtK0s[fNK0s] = pt;
2232  fNK0s++;
2233  }
2234  }
2235  }
2236  }
2237 }
Int_t charge
Bool_t fModulationDEDx
QA pid object.
double Double_t
Definition: External.C:58
AliFlowVZEROResults * fContAllChargesV0Cv3
results
Int_t fNpiPos
pt of K0s in my private selection
Float_t fPhiK0s[1000]
number of K0s in my private selection
void SetDetOR(Int_t idet)
Definition: External.C:236
TProfile * fHResTPCv0C3
also for v3
void SetVarRange(Int_t ivar, Float_t xMin, Float_t xMax)
TProfile * fHResMA3
TProfile for subevent resolution (output)
Float_t fWidthQv3[nCentrBin][2][2]
also for v3
static Float_t fgPsi3tpcMC
static Float_t fgPsi2v0c
static Float_t fgPsi2v0a
Float_t evPlAngV0ACor2
current centrality for the tree
TList * fList2
List for output objects.
static Bool_t fgIsPsiComputed
TProfile * fHResAC2
TProfile for subevent resolution (output)
TProfile * fHResTPCv0A2
...
static Float_t fgPsi2tpcMC
Double_t mass
Float_t fPtK0s[1000]
phi of K0s in my private selection
Float_t GetTOFMismProb() const
Float_t evPlAngV0CCor2
subevent EPs (v2)
centrality
void AddSpecies(const char *name, Int_t nXbin, const Double_t *xbin, Int_t nYbin, const Double_t *ybin)
Int_t fNpiNeg
number of positive pions for K0s selection
virtual void UserCreateOutputObjects()
void ComputeProb(const AliESDtrack *t, Float_t)
virtual void Terminate(Option_t *)
virtual void Analyze(AliAODEvent *aodEvent, Float_t v0Centr)
Float_t evPlAng2
subevent EPs (v2)
void OpenInfoCalbration(Int_t run)
void Fill(Int_t species, Float_t pt, Float_t v2, Float_t x[])
void SetPsiCorrectionDeDx(Float_t psi, Float_t res)
void SetDetResponse(AliESDEvent *esd, Float_t centrality=-1.0, EStartTimeType_t flagStart=AliESDpid::kTOF_T0, Bool_t=kFALSE)
TProfile * fHResMA2
results
Int_t PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD, Int_t specie)
static Float_t fgPsi2tpc
TH2F * fHctauAt1EP
ctau vs DeltaPhi,pt
void Fill(Int_t species, Float_t var1, Float_t var2, Float_t x[])
void SetMC(Bool_t flag)
void QA(UInt_t run, const char *xmlfile="wn.xml", Int_t stage=0, const char *cdb="raw://")
Definition: QA.C:626
Float_t GetExpDeDx(const AliVTrack *t, Int_t iS) const
TH2D * fHdcaPt[nCentrBin][7]
void SetVarName(Int_t ivar, const char *name)
Float_t fWidthQ[nCentrBin][2][2]
and recentering
static Float_t fgPsi3tpc
static Float_t fgPsi2v0cMC
TProfile * fHResv0Cv0A2
TProfile for subevent resolution (output)
TList * fList4
List for output objects.
int Int_t
Definition: External.C:63
AliFlowVZEROResults * fContAllChargesMCC
results
TH2D * fHK0vsLambda
K0s mass vs. pt (standard selection)
TProfile * fHResTPCv0C2
TProfile for subevent resolution (output)
float Float_t
Definition: External.C:68
Bool_t GetCurrentMask(Int_t idet) const
Float_t evPlAng3
subevent EPs (v3)
TTree * fTree
PID class for the Bayesian probabilities.
AliFlowVZEROResults * fContAllChargesMC
results
void SetVarName(Int_t ivar, const char *name)
TH2D * fHK0sMass
EP distribution.
TH2D * fHdcaPtSec[nCentrBin][7]
DCA distribution (for MC primary)
TProfile * fHResMC3
also for v3
TH2D * fHKsPhiEP
Ks phi distribution.
Definition: External.C:228
Int_t fNK0s
primary vertex z coordinate
AliFlowVZEROResults * fContAllChargesMCAv3
results
Double_t nsigma
Int_t fIPiPos[1000]
number of negative pions for K0s selection
void ResetDetOR(Int_t idet)
static Float_t fgPsi2v0aMC
AliFlowVZEROResults * fContAllChargesV0A
subevent EPs (v3)
TList * fList3
List for output objects.
TH2F * fPhiRPv0C
EP distribution vs. centrality (v2)
AliFlowVZEROQA * fQAv3
QA histos (v2)
TProfile * fHResAC3
also for v3
Float_t fMeanQ[nCentrBin][2][2]
loaded by OADB
void AddSpecies(const char *name, Int_t nXbin, const Double_t *bin)
AliFlowVZEROQA * fQA
EP distribution vs. centrality (v3)
TH2F * fPhiRPv0Cv3
EP distribution vs. centrality (v3)
Float_t evPlAngV0CCor3
subevent EPs (v3)
Int_t fNcluster
current run checked to load VZERO calibrations
TProfile * fHResv0Cv0A3
also for v3
AliFlowVZEROResults * fContAllChargesMCCv3
results
Int_t fRun
cenrality bins
Float_t evPlAngV0ACor3
subevent EPs (v2)
Int_t fIPiNeg[1000]
position in the AOD stack of positive pions for K0s
AliFlowVZEROQA * fQA2v3
QA histos (v3)
void SetNewTrackParam(Bool_t flag=kTRUE)
static Float_t fgPsi3v0aMC
void SetVarRange(Int_t ivar, Float_t xMin, Float_t xMax)
AliFlowVZEROResults * fContAllChargesV0C
results
AliFlowVZEROResults * fContAllChargesMCA
also for v3
TProfile * fHResMC2
TProfile for subevent resolution (output)
TProfile * fHResTPCv0A3
TProfile for subevent resolution (output)
Float_t GetDeDx() const
Float_t fCentrality
tree to debug EP (if needed)
static Float_t fgPsi3v0cMC
TH2D * fHK0sMass2
K0s mass vs. pt (private selection)
AliFlowVZEROResults * fContAllChargesV0Av3
results
virtual void UserExec(Option_t *option)
const char Option_t
Definition: External.C:48
TProfile2D * fHctauPtEP
K0s vs lambda mass (in private K0s selection)
Float_t fV0Apol
loaded by OADB
bool Bool_t
Definition: External.C:53
AliFlowVZEROResults * fContQApid
DCA distribution (for MC secondary, not used for data)
TH2F * fPhiRPv0A
also for v3
AliESDtrackCuts * fCutsDaughter
ctau vs. DeltaPhi at 1 GeV/c
AliFlowBayesianPID * fPID
QA histos (v3)
AliESDpid * GetESDpid()
static Float_t fgPsi3v0c
virtual Float_t GetVertex(AliAODEvent *aod) const
AliFlowVZEROQA * fQA2
QA histos (v2)
TH2F * fPhiRPv0Av3
EP distribution vs. centrality (v2)
Float_t fMeanQv3[nCentrBin][2][2]
...
static Float_t fgPsi3v0a
Float_t fV0Cpol
object containing VZERO calibration information
TH2D * fHKsPhi
position in the AOD stack of negative pions for K0s
TProfile * fMultV0
List for output objects.
static const Int_t nCentrBin