]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGPP/pid/AliAnalysisTaskK0sBayes.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
minor fix
[u/mrichter/AliRoot.git] / PWGPP / pid / AliAnalysisTaskK0sBayes.cxx
CommitLineData
a8ad4709 1#include "AliAnalysisTaskK0sBayes.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 "AliESDtrack.h"
12#include "AliCentrality.h"
13#include "AliVHeader.h"
14#include "AliAODVZERO.h"
15#include "TFile.h"
16#include "AliOADBContainer.h"
17#include "TH2F.h"
18#include "TF1.h"
19#include "AliGenHijingEventHeader.h"
20#include "AliMCEvent.h"
21#include "AliAODMCHeader.h"
22#include "AliAODMCParticle.h"
23#include "TChain.h"
24#include "AliESDtrackCuts.h"
25#include "AliESDVertex.h"
26#include "AliEventplane.h"
27#include "AliAnalysisManager.h"
28#include "TRandom.h"
29
30
31Float_t AliAnalysisTaskK0sBayes::fPtKsMin[AliAnalysisTaskK0sBayes::nPtBin] = {0.,0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.};// ptmin bin
32Float_t AliAnalysisTaskK0sBayes::fPtKsMax[AliAnalysisTaskK0sBayes::nPtBin] = {0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4.,5.,6.,10.};// ptmax bin
33
34// STL includes
35//#include <iostream>
36//using namespace std;
37
38ClassImp(AliAnalysisTaskK0sBayes)
39//_____________________________________________________________________________
40AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes():
41 AliAnalysisTaskSE(),
42 fVtxCut(10.0), // cut on |vertex| < fVtxCut
43 fEtaCut(0.8), // cut on |eta| < fEtaCut
44 fMinPt(0.15), // cut on pt > fMinPt
45 fIsMC(kFALSE),
46 fQAsw(kFALSE),
47 fNcluster(70),
48 fFilterBit(1),
49 fList(new TList()),
50 fList2(new TList()),
51 fList3(new TList()),
52 fCentrality(-1),
53 fPsi(0),
54 fPtKs(0.),
55 fPhiKs(0.),
56 fEtaKs(0.),
57 fMassV0(0.),
58 fPtKp(0.),
59 fPhiKp(0.),
60 fEtaKp(0.),
61 fPtKn(0.),
62 fPhiKn(0.),
63 fEtaKn(0.),
64 fPidKp(0),
65 fPidKn(0),
66 fTOFTPCsignal(0),
67 fCombsignal(0),
68 fCutsDaughter(NULL),
69 fPIDCombined(NULL),
70 fContPid(NULL),
71 fContPid2(NULL),
72 fNK0s(0),
73 fNpiPos(0),
74 fNpiNeg(0),
75 fHmismTOF(0),
9c668341 76 fHchannelTOFdistr(0),
77 fTypeCol(2)
a8ad4709 78{
79 // Default constructor (should not be used)
80 fList->SetName("contKsBayes1");
81 fList2->SetName("contKsBayes2");
82 fList3->SetName("contKsBayes3");
83
84 fList->SetOwner(kTRUE);
85 fList2->SetOwner(kTRUE);
86 fList3->SetOwner(kTRUE);
87
88 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
89 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
90
91 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
92 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
f4c4ac9e 93
94 for(Int_t i=0;i < nCentrBin;i++){
95 fElTOF[i] = NULL;
96 fElTPC[i] = NULL;
97 fPiTOF[i] = NULL;
98 fPiTPC[i] = NULL;
99 fKaTOF[i] = NULL;
100 fKaTPC[i] = NULL;
101 fPrTOF[i] = NULL;
102 fPrTPC[i] = NULL;
103 }
104 for(Int_t i=0;i < 4;i++){
105 hMatching[i] = NULL;
106 hTracking[i] = NULL;
107 }
108 for(Int_t i=0;i < 1000;i++){
109 fPhiK0s[i] = 0.0;
110 fPtK0s[i] = 0.0;
111 fIPiPos[i] = 0;
112 fIPiNeg[i] = 0;
113 fIpiP[i] = 0;
114 fIpiN[i] = 0;
115 fMassKs[i] = 0.0;
116 }
a8ad4709 117}
118
119//______________________________________________________________________________
120AliAnalysisTaskK0sBayes::AliAnalysisTaskK0sBayes(const char *name):
121 AliAnalysisTaskSE(name),
122 fVtxCut(10.0), // cut on |vertex| < fVtxCut
123 fEtaCut(0.8), // cut on |eta| < fEtaCut
124 fMinPt(0.15), // cut on pt > fMinPt
125 fIsMC(kFALSE),
126 fQAsw(kFALSE),
127 fNcluster(70),
128 fFilterBit(1),
129 fList(new TList()),
130 fList2(new TList()),
131 fList3(new TList()),
132 fCentrality(-1),
133 fPsi(0),
134 fPtKs(0.),
135 fPhiKs(0.),
136 fEtaKs(0.),
137 fMassV0(0.),
138 fPtKp(0.),
139 fPhiKp(0.),
140 fEtaKp(0.),
141 fPtKn(0.),
142 fPhiKn(0.),
143 fEtaKn(0.),
144 fPidKp(0),
145 fPidKn(0),
146 fTOFTPCsignal(0),
147 fCombsignal(0),
148 fCutsDaughter(NULL),
149 fPIDCombined(NULL),
150 fContPid(NULL),
151 fContPid2(NULL),
152 fNK0s(0),
153 fNpiPos(0),
154 fNpiNeg(0),
155 fHmismTOF(0),
9c668341 156 fHchannelTOFdistr(0),
157 fTypeCol(2)
a8ad4709 158{
159
160 DefineOutput(1, TList::Class());
161 DefineOutput(2, TList::Class());
162 DefineOutput(3, TList::Class());
163 DefineOutput(4, TList::Class());
164
165 // Output slot #1 writes into a TTree
166 fList->SetName("contKsBayes1");
167 fList2->SetName("contKsBayes2");
168 fList3->SetName("contKsBayes3");
169
170 fList->SetOwner(kTRUE);
171 fList2->SetOwner(kTRUE);
172 fList3->SetOwner(kTRUE);
173
174 TFile *fmism = new TFile("$ALICE_ROOT/TOF/data/TOFmismatchDistr.root");
175 fHmismTOF = (TH1F *) fmism->Get("TOFmismDistr");
176
177 TFile *fchDist = new TFile("$ALICE_ROOT/TOF/data/TOFchannelDist.root");
178 fHchannelTOFdistr = (TH1D *) fchDist->Get("hTOFchanDist");
f4c4ac9e 179
180 for(Int_t i=0;i < nCentrBin;i++){
181 fElTOF[i] = NULL;
182 fElTPC[i] = NULL;
183 fPiTOF[i] = NULL;
184 fPiTPC[i] = NULL;
185 fKaTOF[i] = NULL;
186 fKaTPC[i] = NULL;
187 fPrTOF[i] = NULL;
188 fPrTPC[i] = NULL;
189 }
190 for(Int_t i=0;i < 4;i++){
191 hMatching[i] = NULL;
192 hTracking[i] = NULL;
193 }
194 for(Int_t i=0;i < 1000;i++){
195 fPhiK0s[i] = 0.0;
196 fPtK0s[i] = 0.0;
197 fIPiPos[i] = 0;
198 fIPiNeg[i] = 0;
199 fIpiP[i] = 0;
200 fIpiN[i] = 0;
201 fMassKs[i] = 0.0;
202 }
a8ad4709 203}
204//_____________________________________________________________________________
205AliAnalysisTaskK0sBayes::~AliAnalysisTaskK0sBayes()
206{
207
208}
209
210//______________________________________________________________________________
211void AliAnalysisTaskK0sBayes::UserCreateOutputObjects()
212{
213
214 fPIDCombined=new AliPIDCombined;
215 fPIDCombined->SetDefaultTPCPriors();
216 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
217
218 Float_t invMmin = 0.497-0.005*8;
219 Float_t invMmax = 0.497+0.005*8;
220
221 const Int_t nBinPid = 14;
222
223 Int_t binPid[nBinPid] = {1/*ptKs*/,8/*EtaPiP*/,20/*pt+*/,1/*pt-*/,5/*P+*/,1/*P-*/,2/*TOFmatch+*/,1/*TOFmatch-*/,2/*istrue*/,4/*Nsigma+*/,1/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
224
225 Int_t binPid2[nBinPid] = {1/*ptKs*/,8/*EtaPiN*/,1/*pt+*/,20/*pt-*/,1/*P+*/,5/*P-*/,1/*TOFmatch+*/,2/*TOFmatch-*/,2/*istrue*/,1/*Nsigma+*/,4/*Nsigma-*/,1/*DeltaPhi+*/,1/*DeltaPhi-*/,1/*Psi*/};
226
227 fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
228 fContPid->SetTitleX("M_{K^{0}_{s}}");
229 fContPid->SetTitleY("centrality (%)");
230 fContPid->SetVarName(0,"p_{T}^{K^{0}_{s}}");
231 fContPid->SetVarRange(0,0,10);
232 fContPid->SetVarName(1,"#eta^{K^{0}_{s}}");
233 fContPid->SetVarRange(1,-0.8,0.8);
234 fContPid->SetVarName(2,"p_{T}^{Kp}");
235 fContPid->SetVarRange(2,0.3,4.3);
236 fContPid->SetVarName(3,"p_{T}^{Kn}");
237 fContPid->SetVarRange(3,0.3,4.3);
238 fContPid->SetVarName(4,"BayesProb^{Kp}");
239 fContPid->SetVarRange(4,0,1.);
240 fContPid->SetVarName(5,"BayesProb^{Kn}");
241 fContPid->SetVarRange(5,0,1.);
242 fContPid->SetVarName(6,"isTOFmatch^{Kp}");
243 fContPid->SetVarRange(6,-0.5,1.5);
244 fContPid->SetVarName(7,"isTOFmatch^{Kn}");
245 fContPid->SetVarRange(7,-0.5,1.5);
246 fContPid->SetVarName(8,"isKsTrue^{Kn}");
247 fContPid->SetVarRange(8,-0.5,1.5);
248 fContPid->SetVarName(9,"N#sigma^{Kp}");
249 fContPid->SetVarRange(9,1.25,6.25);
250 fContPid->SetVarName(10,"N#sigma^{Kn}");
251 fContPid->SetVarRange(10,1.25,6.25);
252 fContPid->SetVarName(11,"#Delta#phi^{Kp}");
253 fContPid->SetVarRange(11,-TMath::Pi(),TMath::Pi());
254 fContPid->SetVarName(12,"#Delta#phi^{Kn}");
255 fContPid->SetVarRange(12,-TMath::Pi(),TMath::Pi());
256 fContPid->SetVarName(13,"#Psi");
257 fContPid->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
258
259 fContPid2 = new AliPIDperfContainer("contPID2",nBinPid,binPid2);
260 fContPid2->SetTitleX("M_{K^{0}_{s}}");
261 fContPid2->SetTitleY("centrality (%)");
262 fContPid2->SetVarName(0,"p_{T}^{K^{0}_{s}}");
263 fContPid2->SetVarRange(0,0,10);
264 fContPid2->SetVarName(1,"#eta^{K^{0}_{s}}");
265 fContPid2->SetVarRange(1,-0.8,0.8);
266 fContPid2->SetVarName(2,"p_{T}^{Kp}");
267 fContPid2->SetVarRange(2,0.3,4.3);
268 fContPid2->SetVarName(3,"p_{T}^{Kn}");
269 fContPid2->SetVarRange(3,0.3,4.3);
270 fContPid2->SetVarName(4,"BayesProb^{Kp}");
271 fContPid2->SetVarRange(4,0,1.);
272 fContPid2->SetVarName(5,"BayesProb^{Kn}");
273 fContPid2->SetVarRange(5,0,1.);
274 fContPid2->SetVarName(6,"isTOFmatch^{Kp}");
275 fContPid2->SetVarRange(6,-0.5,1.5);
276 fContPid2->SetVarName(7,"isTOFmatch^{Kn}");
277 fContPid2->SetVarRange(7,-0.5,1.5);
278 fContPid2->SetVarName(8,"isKsTrue^{Kn}");
279 fContPid2->SetVarRange(8,-0.5,1.5);
280 fContPid2->SetVarName(9,"N#sigma^{Kp}");
281 fContPid2->SetVarRange(9,1.25,6.25);
282 fContPid2->SetVarName(10,"N#sigma^{Kn}");
283 fContPid2->SetVarRange(10,1.25,6.25);
284 fContPid2->SetVarName(11,"#Delta#phi^{Kp}");
285 fContPid2->SetVarRange(11,-TMath::Pi(),TMath::Pi());
286 fContPid2->SetVarName(12,"#Delta#phi^{Kn}");
287 fContPid2->SetVarRange(12,-TMath::Pi(),TMath::Pi());
288 fContPid2->SetVarName(13,"#Psi");
289 fContPid2->SetVarRange(13,-TMath::Pi()/2,TMath::Pi()/2);
290
291
292
293 const Int_t nDETsignal = 100; // mass
294 Double_t binDETsignal[nDETsignal+1];
295 for(Int_t i=0;i<nDETsignal+1;i++){
296 binDETsignal[i] = invMmin + i*(invMmax - invMmin) / nDETsignal;
297 }
298 const Int_t nDETsignal2 = 10; // centrality
299 Double_t binDETsignal2[nDETsignal2+1];
300 for(Int_t i=0;i<nDETsignal2+1;i++){
301 binDETsignal2[i] = i*100./ nDETsignal2;
302 }
303 fContPid->AddSpecies("K0s",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
304 fContPid2->AddSpecies("K0s2",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
305
306 fList->Add(fContPid);
307 fList->Add(fContPid2);
308
309 hMatching[0] = new TH2F("hMatchAll","TOF matched (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
310 hMatching[1] = new TH2F("hMatchPi","TOF matched (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
311 hMatching[2] = new TH2F("hMatchKa","TOF matched (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
312 hMatching[3] = new TH2F("hMatchPr","TOF matched (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
313
314 hTracking[0] = new TH2F("hTrackingAll","TPC tracks (all);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
315 hTracking[1] = new TH2F("hTrackingPi","TPC tracks (#pi);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
316 hTracking[2] = new TH2F("hTrackingKa","TPC tracks (K);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
317 hTracking[3] = new TH2F("hTrackingPr","TPC tracks (p);p_{T} (GeV/#it{c});centrality (%)",50,0,10,nDETsignal2,0,100);
318
319 fList2->Add(hMatching[0]);
320 fList2->Add(hMatching[1]);
321 fList2->Add(hMatching[2]);
322 fList2->Add(hMatching[3]);
323 fList2->Add(hTracking[0]);
324 fList2->Add(hTracking[1]);
325 fList2->Add(hTracking[2]);
326 fList2->Add(hTracking[3]);
327
328 fTOFTPCsignal = new TH2F("hTOFTPCsignal","TOF-TPC signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma_{TOF};N#sigma_{TPC}",100,-5,5,100,-5,5);
329 fList2->Add(fTOFTPCsignal);
330 fCombsignal = new TH1F("hCombsignal","Combined signal for pions (0.8 < p_{T} < 0.9 GeV/#it{c}, cent. = 0-20);N#sigma",100,0,5);
331 fList2->Add(fCombsignal);
332
333 // QA plots
334 char name[100],title[400];
335 for(Int_t i=0;i < nCentrBin;i++){
336 snprintf(name,100,"hPiTPCQA_%i",i);
337 snprintf(title,400,"TPC signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
338 fPiTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
339 fList3->Add(fPiTPC[i]);
340
341 snprintf(name,100,"hKaTPCQA_%i",i);
342 snprintf(title,400,"TPC signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
343 fKaTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
344 fList3->Add(fKaTPC[i]);
345
346 snprintf(name,100,"hPrTPCQA_%i",i);
347 snprintf(title,400,"TPC signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
348 fPrTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
349 fList3->Add(fPrTPC[i]);
350
351 snprintf(name,100,"hElTPCQA_%i",i);
352 snprintf(title,400,"TPC signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
353 fElTPC[i] = new TH2F(name,title,50,0,4,200,-10,10);
354 fList3->Add(fElTPC[i]);
355
356 snprintf(name,100,"hPiTOFQA_%i",i);
357 snprintf(title,400,"TOF signal for #pi cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
358 fPiTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
359 fList3->Add(fPiTOF[i]);
360
361 snprintf(name,100,"hKaTOFQA_%i",i);
362 snprintf(title,400,"TOF signal for K cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
363 fKaTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
364 fList3->Add(fKaTOF[i]);
365
366 snprintf(name,100,"hPrTOFQA_%i",i);
367 snprintf(title,400,"TOF signal for p cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
368 fPrTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
369 fList3->Add(fPrTOF[i]);
370
371 snprintf(name,100,"hElTOFQA_%i",i);
372 snprintf(title,400,"TOF signal for e cent=%i-%i%c;p_{T} (GeV/#it{c});N#sigma",i*10,(i+1)*10,'%');
373 fElTOF[i] = new TH2F(name,title,50,0,4,200,-10,10);
374 fList3->Add(fElTOF[i]);
375
376 }
377
378 // Post output data.
379 PostData(1, fList);
380 PostData(2, fList2);
381 PostData(3, fList3);
382}
383
384//______________________________________________________________________________
385void AliAnalysisTaskK0sBayes::UserExec(Option_t *)
386{
387 // Main loop
388 // Called for each event
389
390 fOutputAOD = dynamic_cast<AliAODEvent*>(InputEvent());
391 if(!fOutputAOD){
392 Printf("%s:%d AODEvent not found in Input Manager",(char*)__FILE__,__LINE__);
393 this->Dump();
394 return;
395 }
396
397 Float_t zvtx = GetVertex(fOutputAOD);
398
399
400
401 //Get the MC object
402 if(fIsMC){
403 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(fOutputAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
404 if (!mcHeader) {
405 AliError("Could not find MC Header in AOD");
406 return;
407 }
408 }
409
410 if (TMath::Abs(zvtx) < fVtxCut) {
411
412 SelectK0s();
413
414 //Centrality
415 Float_t v0Centr = -10.;
416 Float_t trkCentr = -10.;
417 AliCentrality *centrality = fOutputAOD->GetCentrality();
418 if (centrality){
419 trkCentr = centrality->GetCentralityPercentile("V0M");
420 v0Centr = centrality->GetCentralityPercentile("TRK");
421 }
422
9c668341 423 if((TMath::Abs(v0Centr - trkCentr) < 5.0 || (fTypeCol!=2)) && v0Centr>0){ // consistency cut on centrality selection
a8ad4709 424 fCentrality = v0Centr;
425 Analyze(fOutputAOD); // Do analysis!!!!
426
427 }
428 }
429
430}
431
432//________________________________________________________________________
433void AliAnalysisTaskK0sBayes::Analyze(AliAODEvent* aodEvent)
434{
435
436 Int_t ntrack = aodEvent->GetNumberOfTracks();
437
438 fPtKp=0.,fPhiKp=0.,fEtaKp=0.;
439 fPtKn=0.,fPhiKn=0.,fEtaKn=0.;
440 fPidKp=0,fPidKn=0;
441 fMassV0=-1;
442
443 TLorentzVector ksV;
444
445 Double_t px,py,pz,E;
446
447 Float_t invMmin = 0.497-0.005*8;
448 Float_t invMmax = 0.497+0.005*8;
449
450 Int_t icentr = 8;
451 if(fCentrality < 0) icentr = 8;
452 else if(fCentrality < 10) icentr = 0;
453 else if(fCentrality < 20) icentr = 1;
454 else if(fCentrality < 30) icentr = 2;
455 else if(fCentrality < 40) icentr = 3;
456 else if(fCentrality < 50) icentr = 4;
457 else if(fCentrality < 60) icentr = 5;
458 else if(fCentrality < 70) icentr = 6;
459 else if(fCentrality < 80) icentr = 7;
460
461 Float_t addMismatchForMC = 0.005;
462 if(fCentrality < 50) addMismatchForMC += 0.01;
463 if(fCentrality < 20) addMismatchForMC += 0.02;
464
9c668341 465 if(fTypeCol == 0) addMismatchForMC = 0.005;
466 else if(fTypeCol == 1) addMismatchForMC *= 0.5;
467
a8ad4709 468 fPsi = 0;
469 /* Compute TPC EP */
470 Double_t Qx2 = 0, Qy2 = 0;
471 Double_t Qx3 = 0, Qy3 = 0;
472 for(Int_t iT = 0; iT < ntrack; iT++) {
473 AliAODTrack* aodTrack = aodEvent->GetTrack(iT);
474
475 if (!aodTrack){
476 continue;
477 }
478
479 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
480
481 if ((TMath::Abs(aodTrack->Eta()) > 0.8) || (aodTrack->Pt() < 0.2) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag)
482 continue;
483
484 Double_t b[2] = {-99., -99.};
485 Double_t bCov[3] = {-99., -99., -99.};
486 if (!aodTrack->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 100., b, bCov))
487 continue;
488
489 if ((TMath::Abs(b[0]) > 3.0) || (TMath::Abs(b[1]) > 2.4))
490 continue;
491
492 Qx2 += TMath::Cos(2*aodTrack->Phi());
493 Qy2 += TMath::Sin(2*aodTrack->Phi());
494 Qx3 += TMath::Cos(3*aodTrack->Phi());
495 Qy3 += TMath::Sin(3*aodTrack->Phi());
496
497 }
498
499 fPsi = TMath::ATan2(Qy2, Qx2)/2.;
500
501 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
502 AliInputEventHandler* inputHandler = (AliInputEventHandler*) (man->GetInputEventHandler());
503 AliPIDResponse *PIDResponse=inputHandler->GetPIDResponse();
504
505 PIDResponse->GetTOFResponse().SetTrackParameter(0,0.);
506 PIDResponse->GetTOFResponse().SetTrackParameter(1,0.);
507 PIDResponse->GetTOFResponse().SetTrackParameter(2,0.018);
508 PIDResponse->GetTOFResponse().SetTrackParameter(3,50.0);
509
510 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
511
f4c4ac9e 512 Double_t probP[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
513 Double_t probN[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
a8ad4709 514 Double_t nSigmaTPC,nSigmaTOF=6,nSigmaTPC2,nSigmaTOF2=6,nSigmaComb,nSigmaComb2;
515
516
517 AliAODMCHeader *mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
518 TClonesArray *mcArray = NULL;
519 if (mcHeader)
520 mcArray = (TClonesArray*)aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName());
521
522 Int_t nmc = 0;
523 if(mcArray)
524 nmc = mcArray->GetEntries();
525
526 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
527 AliAODTrack* track = aodEvent->GetTrack(i);
528
529 AliAODMCParticle *mcp = NULL;
530 Int_t pdg = 0;
531
532 if (!track){
533 continue;
534 }
535
536 Int_t tofMatch = (track->GetStatus() & AliVTrack::kTOFout) && (track->GetStatus() & AliVTrack::kTIME);
537
538 Int_t label = -1;
539 if(mcArray){
540 label = track->GetLabel();
541 if(label != -1 && label < nmc){
542 label = TMath::Abs(label);
543 mcp = (AliAODMCParticle*)mcArray->At(label);
544 pdg = TMath::Abs(mcp->GetPdgCode());
545 }
546 else
547 label = -1;
548 }
549 else{
550 /*UInt_t detUsed =*/ fPIDCombined->ComputeProbabilities(track, PIDResponse, probP);
551 }
552
553 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
554 hTracking[0]->Fill(track->Pt(),fCentrality);
555 if(pdg == 211)
556 hTracking[1]->Fill(track->Pt(),fCentrality);
557 else if(pdg == 321)
558 hTracking[2]->Fill(track->Pt(),fCentrality);
559 else if(pdg == 2212)
560 hTracking[3]->Fill(track->Pt(),fCentrality);
561 else if(! mcp){ // Fill matching histo with the prob
562 hTracking[1]->Fill(track->Pt(),fCentrality,probP[2]);
563 hTracking[2]->Fill(track->Pt(),fCentrality,probP[3]);
564 hTracking[3]->Fill(track->Pt(),fCentrality,probP[4]);
565 }
566 }
567
568 if(!tofMatch) continue;
569
570 if(track->TestFilterBit(fFilterBit) && TMath::Abs(track->Eta()) < 0.8 && track->Charge() > 0){
571 hMatching[0]->Fill(track->Pt(),fCentrality);
572 if(pdg == 211)
573 hMatching[1]->Fill(track->Pt(),fCentrality);
574 else if(pdg == 321)
575 hMatching[2]->Fill(track->Pt(),fCentrality);
576 else if(pdg == 2212)
577 hMatching[3]->Fill(track->Pt(),fCentrality);
578 else if(! mcp){ // Fill matching histo with the prob
579 hMatching[1]->Fill(track->Pt(),fCentrality,probP[2]);
580 hMatching[2]->Fill(track->Pt(),fCentrality,probP[3]);
581 hMatching[3]->Fill(track->Pt(),fCentrality,probP[4]);
582 }
583 }
584 }
585
007df888 586// Int_t pdg1 = -1;
587// Int_t pdg2 = -1;
a8ad4709 588
589
590 // start analysis K0s
591 for(Int_t i=0;i < ntrack;i++){ // loop on positive tracks
592 AliAODTrack* KpTrack = aodEvent->GetTrack(i);
593
594 if (!KpTrack){
595 continue;
596 }
597
598 if(!(KpTrack->Charge() > 0 && KpTrack->Pt() > 0.3 && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
599
600 nSigmaComb=5;
e34b28fe 601 nSigmaTOF = 5;
a8ad4709 602 fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
603 fPidKp=0;
604
605 UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
606
607 Double_t oldpP[10];
608 fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
609
610 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
611 fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
612 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
613 fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
614 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
615 fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
616
617 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
618 fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
619
620 if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
621
622 Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
623
9c668341 624 /*
a8ad4709 625 if(mcArray){
626 Int_t labelK = TMath::Abs(KpTrack->GetLabel());
627 AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 628 pdg1 = TMath::Abs(mcp1->GetPdgCode());
a8ad4709 629 }
9c668341 630 */
a8ad4709 631
632 fPidKp = Int_t(probP[2]*100);
633
634 if(tofMatch1){
635 if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
636 // remove this tof hit
637 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
638 detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
639 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
640 fPidKp = Int_t(probP[4]*100);
641 tofMatch1=0;
642 }
643 else{
644 if(probP[2] > probP[3] && probP[2] > probP[4] && probP[2] > probP[0]) fPidKp += 128; // max prob
645
646 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
647 fPrTOF[icentr]->Fill(fPtKp,nSigmaTOF);
648 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
649 fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
650
651 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
652 fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
653
654 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
655 fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
63d1d548 656
a8ad4709 657 if(fIsMC){
658 Float_t mismAdd = addMismatchForMC;
659 if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
660
661 if(gRandom->Rndm() < mismAdd){
662 Float_t etaAbs = TMath::Abs(KpTrack->Eta());
663 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
664 channel = channel % 8736;
665 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
666
667 // generate random time
668 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
669 Double_t times[10];
670 KpTrack->GetIntegratedTimes(times);
671 nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
672 }
673 }
674
675 if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
63d1d548 676 nSigmaTOF = TMath::Abs(nSigmaTOF);
a8ad4709 677
678 if(nSigmaTOF < 2) fPidKp += 256;
679 else if(nSigmaTOF < 3) fPidKp += 512;
680 }
681 }
682
683 if(tofMatch1){
684 nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
685 if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
686 fCombsignal->Fill(nSigmaComb);
687 }
688 } else {
689 nSigmaComb = TMath::Abs(nSigmaTPC);
690 }
691
692 // use sigmaTOF instead of sigmaComb
e34b28fe 693 if(tofMatch1) nSigmaComb = nSigmaTOF;
a8ad4709 694
695 if(nSigmaComb < 2) nSigmaComb = 2;
696 else if(nSigmaComb < 3) nSigmaComb = 3;
697 else if(nSigmaComb < 5) nSigmaComb = 4.99;
698 else nSigmaComb = 6;
699
700 Int_t iks=-1;
701 for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
702 if(i == fIpiP[k]) iks = k;
703 }
704
705 if(fPtKp > 4.299) fPtKp = 4.299;
706
707 if(iks > -1 && fIpiN[iks] > -1){
708 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
709 Int_t j = fIpiN[iks];
710 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
711
712 if (!KnTrack){
713 continue;
714 }
715
716 if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
717
718 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
719 fPidKn=0;
720
721 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
722 Double_t oldpN[10];
723 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
724
725 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
726
727 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
728
729 nSigmaComb2=5;
e34b28fe 730 nSigmaTOF2=5;
a8ad4709 731
732 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 733 /*
a8ad4709 734 if(mcArray){
735 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
736 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 737 pdg2 = TMath::Abs(mcp2->GetPdgCode());
738 }
739 */
a8ad4709 740
741 fPidKn = Int_t(probN[2]*100);
742
743 if(tofMatch2){
744 if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
745 // remove this tof hit
746 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
747 detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
748 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
749 fPidKn = Int_t(probN[4]*100);
750 tofMatch2=0;
751 }
752 else{
753 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
754
755 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
756
757 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
758
759 if(fIsMC){
760 Float_t mismAdd = addMismatchForMC;
761 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
762
763 if(gRandom->Rndm() < mismAdd){
764 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
765 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
766 channel = channel % 8736;
767 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
768
769 // generate random time
770 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
771 Double_t times[10];
772 KnTrack->GetIntegratedTimes(times);
e34b28fe 773 nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
a8ad4709 774 }
775 }
776
777 if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
778
779 if(nSigmaTOF2 < 2) fPidKn += 256;
780 else if(nSigmaTOF2 < 3) fPidKn += 512;
781 }
782 }
783
784 px = KpTrack->Px() + KnTrack->Px();
785 py = KpTrack->Py() + KnTrack->Py();
786 pz = KpTrack->Pz() + KnTrack->Pz();
787 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
788 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
789
790 ksV.SetPxPyPzE(px,py,pz,E);
791 fMassV0 = fMassKs[iks];
792
793 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
794
795
796 fPtKs = ksV.Pt();
797 fEtaKs = ksV.Eta();
798 fPhiKs = ksV.Phi();
799
800 if(tofMatch2){
801 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
802 if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
803 fCombsignal->Fill(nSigmaComb2);
804 }
805 } else {
806 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
807 }
808
809 // use sigmaTOF instead of sigmaComb
e34b28fe 810 if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
a8ad4709 811
812 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
813 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
814 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
815 else nSigmaComb2 = 6;
816
817 Bool_t isTrue = kFALSE;
818
819 if(mcArray){
820 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
821 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
822
823 if(labelP > -1 && labelN > -1){
824 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
825 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
826
827 Int_t mP = partP->GetMother();
828 Int_t mN = partN->GetMother();
829 if(mP == mN && mP > -1){
830 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
831 Int_t pdgM = partM->GetPdgCode();
832 if(pdgM == 310) isTrue = kTRUE;
833 }
834 }
835
836 }
837
838 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
839 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
840
841 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
842 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
843 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
844 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
845
846 if(fPtKn > 4.299) fPtKn = 4.299;
847
848 Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
849
850
851 Int_t ipt = 0;
852 while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
853 ipt++;
854 }
855 ipt--;
856 if(ipt < 0) ipt = 0; // just to be sure
857
858 if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
859 fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
860 xTOfill[1] = KnTrack->Eta();
861 fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
862 }
863
864
865 }
866 } // end analysi K0s
867}
868
869//_____________________________________________________________________________
870Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
871{
872
873 Float_t zvtx = -999;
874
875 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
876 if (!vtxAOD)
877 return zvtx;
878 if(vtxAOD->GetNContributors()>0)
879 zvtx = vtxAOD->GetZ();
880
881 return zvtx;
882}
883//_____________________________________________________________________________
884void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
885{
886 // Terminate loop
887 Printf("Terminate()");
888}
889//=======================================================================
890void AliAnalysisTaskK0sBayes::SelectK0s(){
891 fNK0s=0;
892 fNpiPos=0;
893 fNpiNeg=0;
894
895 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
896 AliAODv0 *myV0;
897 Double_t dMASS=0.0;
898 for (Int_t i=0; i!=nV0s; ++i) {
899 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
900 if(!myV0) continue;
901 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
902 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
903 if(pass) {
904 dMASS = myV0->MassK0Short();
f4c4ac9e 905 Float_t massLambda = myV0->MassLambda();
906 Float_t massAntiLambda = myV0->MassAntiLambda();
907
908 if(TMath::Abs(dMASS-0.497)/0.005 < 8 && TMath::Abs(massLambda-1.115)/0.005 > 8 && TMath::Abs(massAntiLambda-1.115)/0.005 > 8){
a8ad4709 909 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
910 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
911 if(iT->Charge()<0){
912 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
913 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
914 }
915 fPhiK0s[fNK0s] = myV0->Phi();
916 fPtK0s[fNK0s] = myV0->Pt();
917 fIpiP[fNK0s] = FindDaugheterIndex(iT);
918 fIpiN[fNK0s] = FindDaugheterIndex(jT);
919 fMassKs[fNK0s] = dMASS;
920 if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
921 fNK0s++;
922 }
923 }
924 }
925
926 /* My V0 code
927 // fill pion stacks
928 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
929 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
930 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
931
932 if (!aodTrack){
933 continue;
934 }
935
936 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
937
938 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
939 continue;
940 }
941
942 Double_t b[2] = {-99., -99.};
943 Double_t bCov[3] = {-99., -99., -99.};
944 if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
945 continue;
946
947 if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
948
949
950 Int_t charge = aodTrack->Charge();
951 if(charge > 0){
952 fIPiPos[fNpiPos] = iT;
953 fNpiPos++;
954 }
955 else{
956 fIPiNeg[fNpiNeg] = iT;
957 fNpiNeg++;
958 }
959 }
960
961 for(Int_t i=0;i < fNpiPos;i++){
962 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
963 AliESDtrack pipE(pip);
964
965 for(Int_t j=0;j < fNpiNeg;j++){
966 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
967 AliESDtrack pinE(pin);
968
969 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
970
971 Double_t pPos[3];
972 Double_t pNeg[3];
973 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
974 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
975
976 Float_t length = (xp+xn)*0.5;
977
978 Float_t pxs = pPos[0] + pNeg[0];
979 Float_t pys = pPos[1] + pNeg[1];
980 Float_t pzs = pPos[2] + pNeg[2];
981 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);
982
983 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
984 Float_t phi = TMath::ATan2(pys,pxs);
985 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
986
987 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
988
989 if(mindist < 0.4&& length > 0.7 && length < 25){
990 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);
991 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);
992
993 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
994 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
995
996 if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
997 fPhiK0s[fNK0s] = phi;
998 fPtK0s[fNK0s] = pt;
999 fIpiP[fNK0s] =fIPiPos[i] ;
1000 fIpiN[fNK0s] = fIPiNeg[j];
1001 fMassKs[fNK0s] = mass;
1002 fNK0s++;
1003 }
1004 }
1005 }
1006 }
1007 */
1008}
1009
1010//=======================================================================
1011Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1012{
1013 if (myV0->GetOnFlyStatus() ) return 0;
1014
1015 //the following is needed in order to evualuate track-quality
1016 AliAODTrack *iT, *jT;
1017 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1018 Double_t pos[3],cov[6];
1019 vV0s->GetXYZ(pos);
1020 vV0s->GetCovarianceMatrix(cov);
1021 const AliESDVertex vESD(pos,cov,100.,100);
1022
1023 // TESTING CHARGE
1024 int iPos, iNeg;
1025 iT=(AliAODTrack*) myV0->GetDaughter(0);
1026 if(iT->Charge()>0) {
1027 iPos = 0; iNeg = 1;
1028 } else {
1029 iPos = 1; iNeg = 0;
1030 }
1031 // END OF TEST
1032
1033 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1034
1035 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1036
1037 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1038 if(!trkFlag) return 0;
1039 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1040 if(!trkFlag2) return 0;
1041
1042 Double_t pvertex[3];
1043 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1044 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1045 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1046
1047 Double_t dDL=myV0->DecayLengthV0( pvertex );
1048 if(dDL < 0.5 || dDL > 25) return 0;
1049
1050 Double_t dDCA=myV0->DcaV0Daughters();
1051 if(dDCA >0.5) return 0;
1052
1053 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1054 if(dCTP < -1) return 0;
1055
1056// AliESDtrack ieT( iT );
1057// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1058// if(TMath::Abs(dD0P) < 0]) return 0;
1059
1060// AliESDtrack jeT( jT );
1061// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1062// if(TMath::Abs(dD0M) < 0) return 0;
1063
1064// Double_t dD0D0=dD0P*dD0M;
1065// if(dD0D0>0) return 0;
1066
1067// Double_t dETA=myV0->Eta();
1068// if(dETA <-0.8) return 0;
1069// if(dETA >0.8) return 0;
1070
1071// Double_t dQT=myV0->PtArmV0();
1072// if(specie==0) if(dQT<???) return 0;
1073
1074 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1075 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1076
1077 if(specie==1 && dALPHA<0) return 2; // antilambda
1078 return 1; // K0s or lambda
1079}
1080//-------------------------------------------------------------------------
1081Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1082 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1083
1084 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1085 AliAODTrack* track = fOutputAOD->GetTrack(i);
1086 if(track == trk) return i;
1087 }
1088
1089 printf("Daughter for %p not found\n",trk);
1090 return -1;
1091}
1092//-------------------------------------------------------------------------
1093Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1094 if(!fIsMC) return 1; // used only on MC
1095
9c668341 1096 if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1097 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1098
1099 if(!(channel%20)) return 0; // 5% additional loss in MC
1100 }
1101
1102 return 1;
1103}