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