]> 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
warning fix
[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),
76 fHchannelTOFdistr(0)
77{
78 // Default constructor (should not be used)
79 fList->SetName("contLambdaBayes1");
80 fList2->SetName("contLambdaBayes2");
81 fList3->SetName("contLambdaBayes3");
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//______________________________________________________________________________
95AliAnalysisTaskLambdaBayes::AliAnalysisTaskLambdaBayes(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 fPtLambdaC(0.),
110 fPhiLambdaC(0.),
111 fEtaLambda(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 fNLambda(0),
128 fNpPos(0),
129 fNpNeg(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("contLambdaBayes1");
141 fList2->SetName("contLambdaBayes2");
142 fList3->SetName("contLambdaBayes3");
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//_____________________________________________________________________________
155AliAnalysisTaskLambdaBayes::~AliAnalysisTaskLambdaBayes()
156{
157
158}
159
160//______________________________________________________________________________
161void AliAnalysisTaskLambdaBayes::UserCreateOutputObjects()
162{
163
164 fPIDCombined=new AliPIDCombined;
165 fPIDCombined->SetDefaultTPCPriors();
166 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
167
168 Float_t invMmin = 1.115-0.005*8;
169 Float_t invMmax = 1.115+0.005*8;
170
171 const Int_t nBinPid = 14;
172
173 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*/};
174
175 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*/};
176
177 fContPid = new AliPIDperfContainer("contPID",nBinPid,binPid);
178 fContPid->SetTitleX("M_{#Lambda}");
179 fContPid->SetTitleY("centrality (%)");
180 fContPid->SetVarName(0,"p_{T}^{#Lambda}");
181 fContPid->SetVarRange(0,0,10);
182 fContPid->SetVarName(1,"#eta^{#Lambda}");
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,"isLambdaTrue^{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_{#Lambda}");
211 fContPid2->SetTitleY("centrality (%)");
212 fContPid2->SetVarName(0,"p_{T}^{#Lambda}");
213 fContPid2->SetVarRange(0,0,10);
214 fContPid2->SetVarName(1,"#eta^{#Lambda}");
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,"isLambdaTrue^{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("Lambda",nDETsignal,binDETsignal,nDETsignal2,binDETsignal2);
254 fContPid2->AddSpecies("Lambda2",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 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);
279 fList2->Add(fTOFTPCsignal);
280 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);
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 AliAnalysisTaskLambdaBayes::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 SelectLambda();
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 AliAnalysisTaskLambdaBayes::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 = 1.115-0.005*8;
398 Float_t invMmax = 1.115+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 Lambda
537 for(Int_t i=0;i < ntrack;i++){ // loop on proton candidate tracks
538 AliAODTrack* KpTrack = aodEvent->GetTrack(i);
539
540 if (!KpTrack){
541 continue;
542 }
543
544 if(!(KpTrack->Pt() > 0.3 && TMath::Abs(KpTrack->Eta()) < 0.8)) continue;
545
546 Bool_t isLambda = (KpTrack->Charge() > 0);
547
548 nSigmaComb=5;
549 fPtKp=KpTrack->Pt(),fPhiKp=KpTrack->Phi(),fEtaKp=KpTrack->Eta();
550 fPidKp=0;
551
552 UInt_t detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
553
554 Double_t oldpP[10];
555 fPIDCombined->GetPriors(KpTrack, oldpP, PIDResponse, detUsedP);
556
557
558 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kKaon);
559 if(isLambda)fKaTPC[icentr]->Fill(fPtKp,nSigmaTPC);
560 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kElectron);
561 if(isLambda) fElTPC[icentr]->Fill(fPtKp,nSigmaTPC);
562 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kPion);
563 if(isLambda)fPiTPC[icentr]->Fill(fPtKp,nSigmaTPC);
564 nSigmaTPC = PIDResponse->NumberOfSigmasTPC(KpTrack,AliPID::kProton);
565 if(isLambda) fPrTPC[icentr]->Fill(fPtKp,nSigmaTPC);
566
567 if(! (TMath::Abs(nSigmaTPC) < 5)) continue;
568
569 Int_t tofMatch1 = (KpTrack->GetStatus() & AliVTrack::kTOFout) && (KpTrack->GetStatus() & AliVTrack::kTIME);
570
571 if(mcArray){
572 Int_t labelK = TMath::Abs(KpTrack->GetLabel());
573 AliAODMCParticle *mcp1 = (AliAODMCParticle*)mcArray->At(labelK);
007df888 574// pdg1 = TMath::Abs(mcp1->GetPdgCode());
a8ad4709 575 }
576
577 fPidKp = Int_t(probP[4]*100);
578
579 if(tofMatch1){
580 if(!IsChannelValid(TMath::Abs(KpTrack->Eta()))){
581 // remove this tof hit
582 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
583 detUsedP = fPIDCombined->ComputeProbabilities(KpTrack, PIDResponse, probP);
584 fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
585 fPidKp = Int_t(probP[4]*100);
586 tofMatch1=0;
587 }
588 else{
589 if(probP[4] > probP[3] && probP[4] > probP[2] && probP[4] > probP[0]) fPidKp += 128; // max prob
590
591 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kElectron);
592 if(isLambda) fElTOF[icentr]->Fill(fPtKp,nSigmaTOF);
593 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kKaon);
594 if(isLambda) fKaTOF[icentr]->Fill(fPtKp,nSigmaTOF);
595 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kPion);
596 if(isLambda) fPiTOF[icentr]->Fill(fPtKp,nSigmaTOF);
597 nSigmaTOF = PIDResponse->NumberOfSigmasTOF(KpTrack,AliPID::kProton);
598 if(isLambda) fPrTOF[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[4])/PIDResponse->GetTOFResponse().GetExpectedSigma(KpTrack->P(), times[3], AliPID::kProton);
617 }
618 }
619
620 if(fCentrality < 20 && KpTrack->Pt() < 1.3 && KpTrack->Pt() > 1.2)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() < 1.3 && KpTrack->Pt() > 1.2){
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 < fNLambda;k++){ // find the Lambda which contains the positive track
646 if(i == fIpP[k]) iks = k;
647 }
648
649 if(fPtKp > 4.299) fPtKp = 4.299;
650
651 if(iks > -1 && fIpN[iks] > -1){
652 //for(Int_t j=0;j < ntrack;j++){ // loop on negative tracks
653 Int_t j = fIpN[iks];
654 AliAODTrack* KnTrack = aodEvent->GetTrack(j);
655
656 if (!KnTrack){
657 continue;
658 }
659
660 if(!(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(probN[2] > probN[3] && probN[2] > probN[4] && probN[2] > probN[0]) fPidKn += 128; // max prob
687
688 nSigmaTOF2 = PIDResponse->NumberOfSigmasTOF(KnTrack,AliPID::kPion);
689
690 nSigmaTOF2 = TMath::Abs(nSigmaTOF2);
691
692 if(nSigmaTOF2 < 2) fPidKn += 256;
693 else if(nSigmaTOF2 < 3) fPidKn += 512;
694 }
695
696 px = KpTrack->Px() + KnTrack->Px();
697 py = KpTrack->Py() + KnTrack->Py();
698 pz = KpTrack->Pz() + KnTrack->Pz();
699 E = TMath::Sqrt(KpTrack->P()*KpTrack->P() + 0.938e-01*0.938e-01);
700 E += TMath::Sqrt(KnTrack->P()*KnTrack->P()+ 1.39e-01*1.39e-01);
701
702 ksV.SetPxPyPzE(px,py,pz,E);
703 fMassV0 = fMassLambda[iks];
704
705 if(fMassV0 < invMmin || fMassV0 > invMmax) continue;
706
707
708 fPtLambdaC = ksV.Pt();
709 fEtaLambda = ksV.Eta();
710 fPhiLambdaC = ksV.Phi();
711
712 if(tofMatch2){
713 nSigmaComb2 = TMath::Sqrt(nSigmaTOF2*nSigmaTOF2+ nSigmaTPC2*nSigmaTPC2);
714 } else {
715 nSigmaComb2 = TMath::Abs(nSigmaTPC2);
716 }
717
718 // use sigmaTOF instead of sigmaComb
719 //nSigmaComb2 = nSigmaTOF2;
720
721 if(nSigmaComb2 < 2) nSigmaComb2 = 2;
722 else if(nSigmaComb2 < 3) nSigmaComb2 = 3;
723 else if(nSigmaComb2 < 5) nSigmaComb2 = 4.99;
724 else nSigmaComb2 = 6;
725
726 Bool_t isTrue = kFALSE;
727
728 if(mcArray){
729 Int_t labelP = TMath::Abs(KpTrack->GetLabel());
730 Int_t labelN = TMath::Abs(KnTrack->GetLabel());
731
732 if(labelP > -1 && labelN > -1){
733 AliAODMCParticle *partP = (AliAODMCParticle*)mcArray->At(labelP);
734 AliAODMCParticle *partN = (AliAODMCParticle*)mcArray->At(labelN);
735
736 Int_t mP = partP->GetMother();
737 Int_t mN = partN->GetMother();
738 if(mP == mN && mP > -1){
739 AliAODMCParticle *partM = (AliAODMCParticle*)mcArray->At(mP);
740 Int_t pdgM = partM->GetPdgCode();
741 if(TMath::Abs(pdgM) == 3122) isTrue = kTRUE;
742 }
743 }
744
745 }
746
747 Double_t deltaphi1 = KpTrack->Phi() - fPsi;
748 Double_t deltaphi2 = KnTrack->Phi() - fPsi;
749
750 while(deltaphi1 > TMath::Pi()) deltaphi1 -= TMath::Pi()*2;
751 while(deltaphi1 < -TMath::Pi()) deltaphi1 += TMath::Pi()*2;
752 while(deltaphi2 > TMath::Pi()) deltaphi2 -= TMath::Pi()*2;
753 while(deltaphi2 < -TMath::Pi()) deltaphi2 += TMath::Pi()*2;
754
755 if(fPtKn > 4.299) fPtKn = 4.299;
756
757 Float_t xTOfill[] = {fPtLambdaC,KpTrack->Eta(),fPtKp,fPtKn,(fPidKp%128)*0.01,(fPidKn%128)*0.01,tofMatch1,tofMatch2,isTrue,nSigmaComb,nSigmaComb2,deltaphi1,deltaphi2,fPsi};
758 Float_t xTOfill2[] = {fPtLambdaC,KpTrack->Eta(),fPtKn,fPtKp,(fPidKn%128)*0.01,(fPidKp%128)*0.01,tofMatch2,tofMatch1,isTrue,nSigmaComb2,nSigmaComb,deltaphi2,deltaphi1,fPsi};
759
760
761 Int_t ipt = 0;
762 while(fPtLambdaMin[ipt] < fPtLambdaC && ipt < nPtBin){
763 ipt++;
764 }
765 ipt--;
766 if(ipt < 0) ipt = 0; // just to be sure
767
768 if(TMath::Abs(fEtaLambda) < 0.8 && fPtKp > 0.3 && fPtKn > 0.3){
769 if(isLambda) fContPid->Fill(0,fMassV0,fCentrality,xTOfill);
770 else fContPid2->Fill(0,fMassV0,fCentrality,xTOfill2);
771 }
772
773
774 }
775 } // end analysi Lambda
776}
777
778//_____________________________________________________________________________
779Float_t AliAnalysisTaskLambdaBayes::GetVertex(AliAODEvent* aod) const
780{
781
782 Float_t zvtx = -999;
783
784 const AliAODVertex* vtxAOD = aod->GetPrimaryVertex();
785 if (!vtxAOD)
786 return zvtx;
787 if(vtxAOD->GetNContributors()>0)
788 zvtx = vtxAOD->GetZ();
789
790 return zvtx;
791}
792//_____________________________________________________________________________
793void AliAnalysisTaskLambdaBayes::Terminate(Option_t *)
794{
795 // Terminate loop
796 Printf("Terminate()");
797}
798//=======================================================================
799void AliAnalysisTaskLambdaBayes::SelectLambda(){
800 fNLambda=0;
801 fNpPos=0;
802 fNpNeg=0;
803
804 Int_t nV0s = fOutputAOD->GetNumberOfV0s();
805 AliAODv0 *myV0;
806 Double_t dMASS=0.0;
807 for (Int_t i=0; i!=nV0s; ++i) {
808 myV0 = (AliAODv0*) fOutputAOD->GetV0(i);
809 if(!myV0) continue;
810 if(myV0->Pt()<0.1 || TMath::Abs(myV0->Eta()) > 0.8) continue; // skipping low momentum
811 Int_t pass = PassesAODCuts(myV0,fOutputAOD,1); // check for Lambda
812 if(pass) {
813 if(pass==1) dMASS = myV0->MassLambda();
814 else dMASS = myV0->MassAntiLambda();
815
816 if(TMath::Abs(dMASS-1.115)/0.005 < 8){
817 AliAODTrack *iT=(AliAODTrack*) myV0->GetDaughter(0); // positive
818 AliAODTrack *jT=(AliAODTrack*) myV0->GetDaughter(1); // negative
819 if(iT->Charge()<0){
820 iT=(AliAODTrack*) myV0->GetDaughter(1); // positive
821 jT=(AliAODTrack*) myV0->GetDaughter(0); // negative
822 }
823 fPhiLambda[fNLambda] = myV0->Phi();
824 fPtLambda[fNLambda] = myV0->Pt();
825 fIpP[fNLambda] = FindDaugheterIndex(iT);
826 fIpN[fNLambda] = FindDaugheterIndex(jT);
827
828 if(pass==2){
829 Int_t itemp = fIpP[fNLambda];
830 fIpP[fNLambda] = fIpN[fNLambda];
831 fIpN[fNLambda] = itemp;
832 }
833
834 fMassLambda[fNLambda] = dMASS;
835 if(fIpP[fNLambda] > -1 && fIpN[fNLambda] > -1){
836 fNLambda++;
837 }
838
839 }
840 }
841 }
842}
843
844//=======================================================================
845Int_t AliAnalysisTaskLambdaBayes::PassesAODCuts(AliAODv0 *myV0, AliAODEvent *tAOD,Int_t specie)
846{
847 if (myV0->GetOnFlyStatus() ) return 0;
848
849 //the following is needed in order to evualuate track-quality
850 AliAODTrack *iT, *jT;
851 AliAODVertex *vV0s = myV0->GetSecondaryVtx();
852 Double_t pos[3],cov[6];
853 vV0s->GetXYZ(pos);
854 vV0s->GetCovarianceMatrix(cov);
855 const AliESDVertex vESD(pos,cov,100.,100);
856
857 // TESTING CHARGE
858 int iPos, iNeg;
859 iT=(AliAODTrack*) myV0->GetDaughter(0);
860 if(iT->Charge()>0) {
861 iPos = 0; iNeg = 1;
862 } else {
863 iPos = 1; iNeg = 0;
864 }
865 // END OF TEST
866
867 iT=(AliAODTrack*) myV0->GetDaughter(iPos); // positive
868
869 jT=(AliAODTrack*) myV0->GetDaughter(iNeg); // negative
870
871 Bool_t trkFlag = iT->TestFilterBit(fFilterBit);
872 if(!trkFlag) return 0;
873 Bool_t trkFlag2 = jT->TestFilterBit(fFilterBit);
874 if(!trkFlag2) return 0;
875
876 Double_t pvertex[3];
877 pvertex[0]=tAOD->GetPrimaryVertex()->GetX();
878 pvertex[1]=tAOD->GetPrimaryVertex()->GetY();
879 pvertex[2]=tAOD->GetPrimaryVertex()->GetZ();
880
881 Double_t dDL=myV0->DecayLengthV0( pvertex );
882 if(dDL < 0.5 || dDL > 25) return 0;
883
884 Double_t dDCA=myV0->DcaV0Daughters();
885 if(dDCA >0.5) return 0;
886
887 Double_t dCTP=myV0->CosPointingAngle( pvertex );
888 if(dCTP < -1) return 0;
889
890// AliESDtrack ieT( iT );
891// Double_t dD0P=ieT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
892// if(TMath::Abs(dD0P) < 0]) return 0;
893
894// AliESDtrack jeT( jT );
895// Double_t dD0M=jeT.GetD(pvertex[0],pvertex[1],tAOD->GetMagneticField());
896// if(TMath::Abs(dD0M) < 0) return 0;
897
898// Double_t dD0D0=dD0P*dD0M;
899// if(dD0D0>0) return 0;
900
901// Double_t dETA=myV0->Eta();
902// if(dETA <-0.8) return 0;
903// if(dETA >0.8) return 0;
904
905// Double_t dQT=myV0->PtArmV0();
906// if(specie==0) if(dQT<???) return 0;
907
908 Double_t dALPHA=myV0->AlphaV0(); // AlphaV0 -> AODRecoDecat::Alpha -> return 1.-2./(1.+QlProng(0)/QlProng(1));
909 if(myV0->ChargeProng(iPos)<0) dALPHA = -dALPHA; // protects for a change in convention
910
911 if(specie==1 && dALPHA<0) return 2; // antilambda
912 return 1; // Ks or lambda
913}
914//-------------------------------------------------------------------------
915Int_t AliAnalysisTaskLambdaBayes::FindDaugheterIndex(AliAODTrack *trk){
916 Int_t ntrack = fOutputAOD->GetNumberOfTracks();
917
918 for(Int_t i=0;i < ntrack;i++){ // loop on tracks
919 AliAODTrack* track = fOutputAOD->GetTrack(i);
920 if(track == trk) return i;
921 }
922
923 printf("Daughter for %p not found\n",trk);
924 return -1;
925}
926//-------------------------------------------------------------------------
927Int_t AliAnalysisTaskLambdaBayes::IsChannelValid(Float_t etaAbs){
928 if(!fIsMC) return 1; // used only on MC
929
930 Int_t run = fOutputAOD->GetRunNumber();
931 if( (run>=136851&&run<=139846) || (run>=165772 && run<=170718) ){ // LHC10h or LHC11h because of TOF matching window at 3 cm
932 Int_t channel = Int_t(4334.09 - 4758.36 * etaAbs -1989.71 * etaAbs*etaAbs + 1957.62*etaAbs*etaAbs*etaAbs);
933
934 if(!(channel%20)) return 0; // 5% additional loss in MC
935 }
936
937 return 1;
938}