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