]> 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
update variable name (G.Luparello)
[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);
656
657 nSigmaTOF = TMath::Abs(nSigmaTOF);
658
659 if(fIsMC){
660 Float_t mismAdd = addMismatchForMC;
661 if(KpTrack->Pt() < 1) mismAdd = addMismatchForMC/KpTrack->Pt();
662
663 if(gRandom->Rndm() < mismAdd){
664 Float_t etaAbs = TMath::Abs(KpTrack->Eta());
665 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
666 channel = channel % 8736;
667 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
668
669 // generate random time
670 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
671 Double_t times[10];
672 KpTrack->GetIntegratedTimes(times);
673 nSigmaTOF = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kPion);
674 }
675 }
676
677 if(fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8)fTOFTPCsignal->Fill(nSigmaTOF,nSigmaTPC);
678
679 if(nSigmaTOF < 2) fPidKp += 256;
680 else if(nSigmaTOF < 3) fPidKp += 512;
681 }
682 }
683
684 if(tofMatch1){
685 nSigmaComb = TMath::Sqrt(nSigmaTOF*nSigmaTOF + nSigmaTPC*nSigmaTPC);
686 if(nSigmaTOF < 5 && fCentrality < 20 && KpTrack->Pt() < 0.9 && KpTrack->Pt() > 0.8){
687 fCombsignal->Fill(nSigmaComb);
688 }
689 } else {
690 nSigmaComb = TMath::Abs(nSigmaTPC);
691 }
692
693 // use sigmaTOF instead of sigmaComb
e34b28fe 694 if(tofMatch1) nSigmaComb = nSigmaTOF;
a8ad4709 695
696 if(nSigmaComb < 2) nSigmaComb = 2;
697 else if(nSigmaComb < 3) nSigmaComb = 3;
698 else if(nSigmaComb < 5) nSigmaComb = 4.99;
699 else nSigmaComb = 6;
700
701 Int_t iks=-1;
702 for(Int_t k=0;k < fNK0s;k++){ // find the K0s which contains the positive track
703 if(i == fIpiP[k]) iks = k;
704 }
705
706 if(fPtKp > 4.299) fPtKp = 4.299;
707
708 if(iks > -1 && fIpiN[iks] > -1){
709 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
710 Int_t j = fIpiN[iks];
711 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
712
713 if (!KnTrack){
714 continue;
715 }
716
717 if(!(KnTrack->Charge() < 0 && KnTrack->Pt() > 0.3 && TMath::Abs(KnTrack->Eta()) < 0.8)) continue;
718
719 fPtKn=KnTrack->Pt(),fPhiKn=KnTrack->Phi(),fEtaKn=KnTrack->Eta();
720 fPidKn=0;
721
722 UInt_t detUsedN = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
723 Double_t oldpN[10];
724 fPIDCombined->GetPriors(KnTrack, oldpN, PIDResponse, detUsedN);
725
726 nSigmaTPC2 = PIDResponse->NumberOfSigmasTPC(KnTrack,AliPID::kPion);
727
728 if(! (TMath::Abs(nSigmaTPC2) < 5)) continue;
729
730 nSigmaComb2=5;
e34b28fe 731 nSigmaTOF2=5;
a8ad4709 732
733 Int_t tofMatch2 = (KnTrack->GetStatus() & AliVTrack::kTOFout) && (KnTrack->GetStatus() & AliVTrack::kTIME);
9c668341 734 /*
a8ad4709 735 if(mcArray){
736 Int_t labelK = TMath::Abs(KnTrack->GetLabel());
737 AliAODMCParticle *mcp2 = (AliAODMCParticle*)mcArray->At(labelK);
9c668341 738 pdg2 = TMath::Abs(mcp2->GetPdgCode());
739 }
740 */
a8ad4709 741
742 fPidKn = Int_t(probN[2]*100);
743
744 if(tofMatch2){
745 if(!IsChannelValid(TMath::Abs(KnTrack->Eta()))){
746 // remove this tof hit
747 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
748 detUsedP = fPIDCombined->ComputeProbabilities(KnTrack, PIDResponse, probN);
749 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
750 fPidKn = Int_t(probN[4]*100);
751 tofMatch2=0;
752 }
753 else{
754 if(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
755
756 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
757
758 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
759
760 if(fIsMC){
761 Float_t mismAdd = addMismatchForMC;
762 if(KnTrack->Pt() < 1) mismAdd = addMismatchForMC/KnTrack->Pt();
763
764 if(gRandom->Rndm() < mismAdd){
765 Float_t etaAbs = TMath::Abs(KnTrack->Eta());
766 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
767 channel = channel % 8736;
768 Float_t distIP = fHchannelTOFdistr->GetBinContent(channel);
769
770 // generate random time
771 Float_t timeRandom = fHmismTOF->GetRandom() + distIP*3.35655419905265973e+00;
772 Double_t times[10];
773 KnTrack->GetIntegratedTimes(times);
e34b28fe 774 nSigmaTOF2 = TMath::Abs(timeRandom - times[2])/PIDResponse->GetTOFResponse().GetExpectedSigma(KnTrack->P(), times[3], AliPID::kPion);
a8ad4709 775 }
776 }
777
778 if(fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1) fTOFTPCsignal->Fill(nSigmaTOF2,nSigmaTPC2);
779
780 if(nSigmaTOF2 < 2) fPidKn += 256;
781 else if(nSigmaTOF2 < 3) fPidKn += 512;
782 }
783 }
784
785 px = KpTrack->Px() + KnTrack->Px();
786 py = KpTrack->Py() + KnTrack->Py();
787 pz = KpTrack->Pz() + KnTrack->Pz();
788 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 1.39e-01*1.39e-01);
789 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
790
791 ksV.SetPxPyPzE(px,py,pz,E);
792 fMassV0 = fMassKs[iks];
793
794 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
795
796
797 fPtKs = ksV.Pt();
798 fEtaKs = ksV.Eta();
799 fPhiKs = ksV.Phi();
800
801 if(tofMatch2){
802 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
803 if(nSigmaTOF2 < 5 && fCentrality < 20 && KnTrack->Pt() < 1.2 && KnTrack->Pt() > 1){
804 fCombsignal->Fill(nSigmaComb2);
805 }
806 } else {
807 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
808 }
809
810 // use sigmaTOF instead of sigmaComb
e34b28fe 811 if(tofMatch2) nSigmaComb2 = nSigmaTOF2;
a8ad4709 812
813 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
814 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
815 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
816 else nSigmaComb2 = 6;
817
818 Bool_t isTrue = kFALSE;
819
820 if(mcArray){
821 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
822 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
823
824 if(labelP > -1 && labelN > -1){
825 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
826 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
827
828 Int_t mP = partP->GetMother();
829 Int_t mN = partN->GetMother();
830 if(mP == mN && mP > -1){
831 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
832 Int_t pdgM = partM->GetPdgCode();
833 if(pdgM == 310) isTrue = kTRUE;
834 }
835 }
836
837 }
838
839 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
840 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
841
842 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
843 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
844 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
845 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
846
847 if(fPtKn > 4.299) fPtKn = 4.299;
848
849 Float_t xTOfill[] = {fPtKs,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
850
851
852 Int_t ipt = 0;
853 while(fPtKsMin[ipt] < fPtKs && ipt < nPtBin){
854 ipt++;
855 }
856 ipt--;
857 if(ipt < 0) ipt = 0; // just to be sure
858
859 if(TMath::Abs(fEtaKs) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
860 fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
861 xTOfill[1] = KnTrack->Eta();
862 fContPid2->Fill(0,fMassV0,fCentrality,xTOfill);
863 }
864
865
866 }
867 } // end analysi K0s
868}
869
870//_____________________________________________________________________________
871Float_t AliAnalysisTaskK0sBayes::GetVertex(AliAODEvent* aod) const
872{
873
874 Float_t zvtx = -999;
875
876 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
877 if (!vtxAOD)
878 return zvtx;
879 if(vtxAOD->GetNContributors()>0)
880 zvtx = vtxAOD->GetZ();
881
882 return zvtx;
883}
884//_____________________________________________________________________________
885void AliAnalysisTaskK0sBayes::Terminate(Option_t *)
886{
887 // Terminate loop
888 Printf("Terminate()");
889}
890//=======================================================================
891void AliAnalysisTaskK0sBayes::SelectK0s(){
892 fNK0s=0;
893 fNpiPos=0;
894 fNpiNeg=0;
895
896 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
897 AliAODv0 *myV0;
898 Double_t dMASS=0.0;
899 for (Int_t i=0; i!=nV0s; ++i) {
900 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
901 if(!myV0) continue;
902 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
903 Int_t pass = PassesAODCuts(myV0,fOutputAOD,0); // check for K0s
904 if(pass) {
905 dMASS = myV0->MassK0Short();
f4c4ac9e 906 Float_t massLambda = myV0->MassLambda();
907 Float_t massAntiLambda = myV0->MassAntiLambda();
908
909 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 910 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
911 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
912 if(iT->Charge()<0){
913 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
914 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
915 }
916 fPhiK0s[fNK0s] = myV0->Phi();
917 fPtK0s[fNK0s] = myV0->Pt();
918 fIpiP[fNK0s] = FindDaugheterIndex(iT);
919 fIpiN[fNK0s] = FindDaugheterIndex(jT);
920 fMassKs[fNK0s] = dMASS;
921 if(fIpiP[fNK0s] > -1 && fIpiN[fNK0s] > -1)
922 fNK0s++;
923 }
924 }
925 }
926
927 /* My V0 code
928 // fill pion stacks
929 Int_t nAODTracks = fOutputAOD->GetNumberOfTracks();
930 for(Int_t iT = 0; iT < nAODTracks; iT++) { // loop on the tracks
931 AliAODTrack* aodTrack = fOutputAOD->GetTrack(iT);
932
933 if (!aodTrack){
934 continue;
935 }
936
937 Bool_t trkFlag = aodTrack->TestFilterBit(fFilterBit);
938
939 if ((TMath::Abs(aodTrack->Eta()) > fEtaCut) || (aodTrack->Pt() < fMinPt) || (aodTrack->GetTPCNcls() < fNcluster) || !trkFlag){
940 continue;
941 }
942
943 Double_t b[2] = {-99., -99.};
944 Double_t bCov[3] = {-99., -99., -99.};
945 if (!aodTrack->PropagateToDCA(fOutputAOD->GetPrimaryVertex(), fOutputAOD->GetMagneticField(), 100., b, bCov))
946 continue;
947
948 if(TMath::Abs(b[0]) < 0.3/aodTrack->Pt()) continue;
949
950
951 Int_t charge = aodTrack->Charge();
952 if(charge > 0){
953 fIPiPos[fNpiPos] = iT;
954 fNpiPos++;
955 }
956 else{
957 fIPiNeg[fNpiNeg] = iT;
958 fNpiNeg++;
959 }
960 }
961
962 for(Int_t i=0;i < fNpiPos;i++){
963 AliAODTrack *pip = fOutputAOD->GetTrack(fIPiPos[i]);
964 AliESDtrack pipE(pip);
965
966 for(Int_t j=0;j < fNpiNeg;j++){
967 AliAODTrack *pin = fOutputAOD->GetTrack(fIPiNeg[j]);
968 AliESDtrack pinE(pin);
969
970 Double_t xn, xp, mindist=pinE.GetDCA(&pipE,fOutputAOD->GetMagneticField(),xn,xp);
971
972 Double_t pPos[3];
973 Double_t pNeg[3];
974 pipE.GetPxPyPzAt(xp,fOutputAOD->GetMagneticField(),pPos);
975 pinE.GetPxPyPzAt(xn,fOutputAOD->GetMagneticField(),pNeg);
976
977 Float_t length = (xp+xn)*0.5;
978
979 Float_t pxs = pPos[0] + pNeg[0];
980 Float_t pys = pPos[1] + pNeg[1];
981 Float_t pzs = pPos[2] + pNeg[2];
982 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);
983
984 Float_t pt = TMath::Sqrt(pxs*pxs + pys*pys);
985 Float_t phi = TMath::ATan2(pys,pxs);
986 Float_t mass = TMath::Sqrt(es*es - pt*pt - pzs*pzs);
987
988 // if(length > 1) printf("length = %f - distance = %f - mass= %f\n",length,mindist,mass);
989
990 if(mindist < 0.4&& length > 0.7 && length < 25){
991 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);
992 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);
993
994 Float_t massaL = TMath::Sqrt(esL*esL - pt*pt - pzs*pzs);
995 Float_t massaAL = TMath::Sqrt(esAL*esAL - pt*pt - pzs*pzs);
996
997 if(TMath::Abs(mass-0.497)/0.005 < 8 && massaL > 1.15 && massaAL > 1.15){
998 fPhiK0s[fNK0s] = phi;
999 fPtK0s[fNK0s] = pt;
1000 fIpiP[fNK0s] =fIPiPos[i] ;
1001 fIpiN[fNK0s] = fIPiNeg[j];
1002 fMassKs[fNK0s] = mass;
1003 fNK0s++;
1004 }
1005 }
1006 }
1007 }
1008 */
1009}
1010
1011//=======================================================================
1012Int_t AliAnalysisTaskK0sBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
1013{
1014 if (myV0->GetOnFlyStatus() ) return 0;
1015
1016 //the following is needed in order to evualuate track-quality
1017 AliAODTrack *iT, *jT;
1018 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
1019 Double_t pos[3],cov[6];
1020 vV0s->GetXYZ(pos);
1021 vV0s->GetCovarianceMatrix(cov);
1022 const AliESDVertex vESD(pos,cov,100.,100);
1023
1024 // TESTING CHARGE
1025 int iPos, iNeg;
1026 iT=(AliAODTrack*) myV0->GetDaughter(0);
1027 if(iT->Charge()>0) {
1028 iPos = 0; iNeg = 1;
1029 } else {
1030 iPos = 1; iNeg = 0;
1031 }
1032 // END OF TEST
1033
1034 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
1035
1036 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
1037
1038 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
1039 if(!trkFlag) return 0;
1040 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
1041 if(!trkFlag2) return 0;
1042
1043 Double_t pvertex[3];
1044 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
1045 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
1046 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
1047
1048 Double_t dDL=myV0->DecayLengthV0( pvertex );
1049 if(dDL < 0.5 || dDL > 25) return 0;
1050
1051 Double_t dDCA=myV0->DcaV0Daughters();
1052 if(dDCA >0.5) return 0;
1053
1054 Double_t dCTP=myV0->CosPointingAngle( pvertex );
1055 if(dCTP < -1) return 0;
1056
1057// AliESDtrack ieT( iT );
1058// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1059// if(TMath::Abs(dD0P) < 0]) return 0;
1060
1061// AliESDtrack jeT( jT );
1062// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
1063// if(TMath::Abs(dD0M) < 0) return 0;
1064
1065// Double_t dD0D0=dD0P*dD0M;
1066// if(dD0D0>0) return 0;
1067
1068// Double_t dETA=myV0->Eta();
1069// if(dETA <-0.8) return 0;
1070// if(dETA >0.8) return 0;
1071
1072// Double_t dQT=myV0->PtArmV0();
1073// if(specie==0) if(dQT<???) return 0;
1074
1075 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
1076 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
1077
1078 if(specie==1 && dALPHA<0) return 2; // antilambda
1079 return 1; // K0s or lambda
1080}
1081//-------------------------------------------------------------------------
1082Int_t AliAnalysisTaskK0sBayes::FindDaugheterIndex(AliAODTrack *trk){
1083 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
1084
1085 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
1086 AliAODTrack* track = fOutputAOD->GetTrack(i);
1087 if(track == trk) return i;
1088 }
1089
1090 printf("Daughter for %p not found\n",trk);
1091 return -1;
1092}
1093//-------------------------------------------------------------------------
1094Int_t AliAnalysisTaskK0sBayes::IsChannelValid(Float_t etaAbs){
1095 if(!fIsMC) return 1; // used only on MC
1096
9c668341 1097 if(fTypeCol == 2){ // LHC10h or LHC11h because of TOF matching window at 3 cm
a8ad4709 1098 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
1099
1100 if(!(channel%20)) return 0; // 5% additional loss in MC
1101 }
1102
1103 return 1;
1104}