]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/STRANGENESS/LambdaK0PbPb/AliAnalysisTaskCTauPbPb.cxx
Support for correlated error calculation in N-dimensional unfolding
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / AliAnalysisTaskCTauPbPb.cxx
CommitLineData
e26966dc 1#include <TCanvas.h>
2#include <TTree.h>
3#include <TFile.h>
4#include <TH1F.h>
5#include <TH2F.h>
6#include <TH3F.h>
7#include <TPDGCode.h>
8#include <TDatabasePDG.h>
9#include <TParticlePDG.h>
10#include <TParticle.h>
11#include <TROOT.h>
12
13#include "AliESDEvent.h"
14#include "AliESDv0.h"
15#include "AliESDcascade.h"
16
17#include "AliCentrality.h"
18
19#include "AliMCEvent.h"
20#include "AliStack.h"
21
22#include "AliPID.h"
23#include "AliPIDResponse.h"
24
25#include "AliInputEventHandler.h"
26#include "AliAnalysisManager.h"
27
28#include "AliAnalysisTaskCTauPbPb.h"
29
30extern TROOT *gROOT;
31
32ClassImp(AliAnalysisTaskCTauPbPb)
33
25aeb948 34static Int_t nbins=100; // number of bins
e26966dc 35static Double_t lMin=0.0, lMax=100.;
36static Double_t pMin=0.0, pMax=10.;
37static Double_t yMax=0.5;
38
39
40//
41// This is a little task for checking the c*tau of the strange particles
42//
43
44AliAnalysisTaskCTauPbPb::AliAnalysisTaskCTauPbPb(const char *name) :
45AliAnalysisTaskSE(name),
46fIsMC(kFALSE),
47fCMin(0.),
48fCMax(90.),
49fOutput(0),
50fMult(0),
51fdEdx(0),
52fdEdxPid(0),
53
54fK0sM(0),
55fK0sSi(0),
56fK0sMC(0),
57fK0sAs(0),
58
59fLambdaM(0),
60fLambdaSi(0),
61fLambdaMC(0),
62fLambdaAs(0),
63
0b384e23 64fLambdaBarM(0),
65fLambdaBarSi(0),
66fLambdaBarMC(0),
67fLambdaBarAs(0),
68
e26966dc 69fCPA(0),
70fDCA(0),
71
72fLambdaEff(0),
73fLambdaPt(0),
74
75fLambdaFromXi(0),
76fXiM(0),
0b384e23 77fXiSiP(0),
78
79fLambdaBarFromXiBar(0),
80fXiBarM(0),
81fXiBarSiP(0)
e26966dc 82{
83 // Constructor. Initialization of pointers
84 DefineOutput(1, TList::Class());
85}
86
87void AliAnalysisTaskCTauPbPb::UserCreateOutputObjects()
88{
89 fOutput = new TList();
90 fOutput->SetOwner();
91
92
93 fMult=new TH1F("fMult","Multiplicity",1100,0.,3300);
94 fMult->GetXaxis()->SetTitle("N tracks");
95 fOutput->Add(fMult);
96
97 fdEdx=new TH2F("fdEdx","dE/dx",50,0.2,3,50,0.,6.);
98 fOutput->Add(fdEdx);
99
100 fdEdxPid=new TH2F("fdEdxPid","dE/dx with PID",50,0.2,3,50,0.,6.);
101 fOutput->Add(fdEdxPid);
102
103 fK0sM =
104 new TH2F("fK0sM", "Mass for K^{0}_{s}", nbins/2, 0.448, 0.548, 10,pMin,pMax);
105 fK0sM->GetXaxis()->SetTitle("Mass [GeV/c]");
106 fOutput->Add(fK0sM);
107
108 fK0sSi =
109 new TH2F("fK0sSi","L_{T} vs p_{T} for K^{0}_{s}, side-band subtracted",
110 nbins,pMin,pMax,nbins,lMin,lMax);
111 fK0sSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
112 fK0sSi->GetYaxis()->SetTitle("L_{T} [cm]");
113 fOutput->Add(fK0sSi);
114
115 fK0sMC =
116 new TH2F("fK0sMC","L_{T} vs p_{T} for K^{0}_{s}, from MC stack",
117 nbins,pMin,pMax,nbins,lMin,lMax);
118 fK0sMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
119 fK0sMC->GetYaxis()->SetTitle("L_{T} [cm]");
120 fOutput->Add(fK0sMC);
121
122 fK0sAs =
123 new TH2F("fK0sAs", "L_{T} vs p_{T} for K^{0}_{s}, associated",
124 nbins,pMin,pMax,nbins,lMin,lMax);
125 fK0sAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
126 fK0sAs->GetYaxis()->SetTitle("L_{T} [cm]");
127 fOutput->Add(fK0sAs);
128
129 //----------------------
130
131 fLambdaM =
132 new TH2F("fLambdaM","Mass for \\Lambda", nbins, 1.065, 1.165,nbins,pMin,pMax);
e26966dc 133 fLambdaM->GetXaxis()->SetTitle("Mass [GeV/c]");
134 fOutput->Add(fLambdaM);
135
136 fLambdaSi =
137 new TH2F("fLambdaSi","L_{T} vs p_{T} for \\Lambda, side-band subtructed",
138 nbins,pMin,pMax,nbins,lMin,lMax);
139 fLambdaSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
140 fLambdaSi->GetYaxis()->SetTitle("L_{T} [cm]");
141 fOutput->Add(fLambdaSi);
142
143 fLambdaMC =
144 new TH2F("fLambdaMC","c\\tau for \\Lambda, from MC stack",
145 nbins,pMin,pMax,nbins,lMin,lMax);
146 fLambdaMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
147 fLambdaMC->GetYaxis()->SetTitle("L_{T} [cm]");
148 fOutput->Add(fLambdaMC);
149
150 fLambdaAs =
151 new TH2F("fLambdaAs","c\\tau for \\Lambda, associated",
152 nbins,pMin,pMax,nbins,lMin,lMax);
153 fLambdaAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
154 fLambdaAs->GetYaxis()->SetTitle("L_{T} [cm]");
155 fOutput->Add(fLambdaAs);
156
0b384e23 157
158 fLambdaBarM =
159 new TH2F("fLambdaBarM","Mass for anti-\\Lambda", nbins, 1.065, 1.165,nbins,pMin,pMax);
160 fLambdaBarM->GetXaxis()->SetTitle("Mass [GeV/c]");
161 fOutput->Add(fLambdaBarM);
162
163 fLambdaBarSi =
164 new TH2F("fLambdaBarSi","L_{T} vs p_{T} for anti-\\Lambda, side-band subtructed",
165 nbins,pMin,pMax,nbins,lMin,lMax);
166 fLambdaBarSi->GetXaxis()->SetTitle("p_{T} [GeV/c]");
167 fLambdaBarSi->GetYaxis()->SetTitle("L_{T} [cm]");
168 fOutput->Add(fLambdaBarSi);
169
170 fLambdaBarMC =
171 new TH2F("fLambdaBarMC","c\\tau for anti-\\Lambda, from MC stack",
172 nbins,pMin,pMax,nbins,lMin,lMax);
173 fLambdaBarMC->GetXaxis()->SetTitle("p_{T} [GeV/c]");
174 fLambdaBarMC->GetYaxis()->SetTitle("L_{T} [cm]");
175 fOutput->Add(fLambdaBarMC);
176
177 fLambdaBarAs =
178 new TH2F("fLambdaBarAs","c\\tau for anti-\\Lambda, associated",
179 nbins,pMin,pMax,nbins,lMin,lMax);
180 fLambdaBarAs->GetXaxis()->SetTitle("p_{T} [GeV/c]");
181 fLambdaBarAs->GetYaxis()->SetTitle("L_{T} [cm]");
182 fOutput->Add(fLambdaBarAs);
183
184
185
e26966dc 186 fCPA=new TH1F("fCPA","Cosine of the pointing angle",30,0.9978,1.);
187 fOutput->Add(fCPA);
188 fDCA=new TH1F("fDCA","DCA between the daughters",30,0.,1.1);
189 fOutput->Add(fDCA);
190
191 fLambdaEff=fLambdaAs->ProjectionX();
192 fLambdaEff->SetName("fLambdaEff");
193 fLambdaEff->SetTitle("Efficiency for #Lambda");
194 fOutput->Add(fLambdaEff);
195
196 fLambdaPt=fLambdaAs->ProjectionX();
197 fLambdaPt->SetName("fLambdaPt");
198 fLambdaPt->SetTitle("Raw #Lambda pT spectrum");
199 fOutput->Add(fLambdaPt);
200
201 //----------------------
202
0b384e23 203 fLambdaFromXi=new TH3F("fLambdaFromXi","L_{T} vs p_{T} vs p_{T} of \\Xi for \\Lambda from \\Xi",
e26966dc 204 nbins,pMin,pMax,nbins,lMin,lMax,33,pMin,pMax+2);
205 fOutput->Add(fLambdaFromXi);
206
207 fXiM =
208 new TH2F("fXiM", "\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
209 fOutput->Add(fXiM);
210
211 fXiSiP = new TH1F("fXiSiP", "Pt for \\Xi, side-band subracted",
212 33,pMin,pMax+2);
213 fOutput->Add(fXiSiP);
214
215
0b384e23 216 fLambdaBarFromXiBar=new TH3F("fLambdaBarFromXiBar","L_{T} vs p_{T} vs p_{T} of anti-\\Xi for anti-\\Lambda from anti-\\Xi",
217 nbins,pMin,pMax,nbins,lMin,lMax,33,pMin,pMax+2);
218 fOutput->Add(fLambdaBarFromXiBar);
219
220 fXiBarM =
221 new TH2F("fXiBarM", "anti-\\Xi mass distribution", 50, 1.271, 1.371,12,pMin,pMax+2);
222 fOutput->Add(fXiBarM);
223
224 fXiBarSiP = new TH1F("fXiBarSiP", "Pt for anti-\\Xi, side-band subracted",
225 33,pMin,pMax+2);
226 fOutput->Add(fXiBarSiP);
227
228
e26966dc 229 PostData(1, fOutput);
230}
231
232static Bool_t AcceptTrack(const AliESDtrack *t) {
233 if (!t->IsOn(AliESDtrack::kTPCrefit)) return kFALSE;
234 if (t->GetKinkIndex(0)>0) return kFALSE;
235
236 Float_t nCrossedRowsTPC = t->GetTPCClusterInfo(2,1);
237 if (nCrossedRowsTPC < 70) return kFALSE;
238 Int_t findable=t->GetTPCNclsF();
239 if (findable <= 0) return kFALSE;
240 if (nCrossedRowsTPC/findable < 0.8) return kFALSE;
241
ff909752 242 if (TMath::Abs(t->Eta()) > 0.8) return kFALSE;
243
e26966dc 244 return kTRUE;
245}
246
247static Bool_t AcceptV0(const AliESDv0 *v0, const AliESDEvent *esd) {
248
249 if (v0->GetOnFlyStatus()) return kFALSE;
250
251 if (v0->Pt() < pMin) return kFALSE;
252
253 Int_t nidx=TMath::Abs(v0->GetNindex());
254 AliESDtrack *ntrack=esd->GetTrack(nidx);
255 if (!AcceptTrack(ntrack)) return kFALSE;
256
257 Int_t pidx=TMath::Abs(v0->GetPindex());
258 AliESDtrack *ptrack=esd->GetTrack(pidx);
259 if (!AcceptTrack(ptrack)) return kFALSE;
260
261 Float_t xy,z0;
262 ntrack->GetImpactParameters(xy,z0);
263 if (TMath::Abs(xy)<0.1) return kFALSE;
264 ptrack->GetImpactParameters(xy,z0);
265 if (TMath::Abs(xy)<0.1) return kFALSE;
266
267 Double_t dca=v0->GetDcaV0Daughters();
268 if (dca>1.0) return kFALSE;
269 //if (dca>0.7) return kFALSE;
270 //if (dca>0.4) return kFALSE;
271
272 Double_t cpa=v0->GetV0CosineOfPointingAngle();
273 if (cpa<0.998) return kFALSE;
274 //if (cpa<0.99875) return kFALSE;
275 //if (cpa<0.9995) return kFALSE;
276
277 Double_t xx,yy,zz; v0->GetXYZ(xx,yy,zz);
278 Double_t r2=xx*xx + yy*yy;
279 if (r2<0.9*0.9) return kFALSE;
280 if (r2>100*100) return kFALSE;
281
282 return kTRUE;
283}
284
285static Bool_t AcceptCascade(const AliESDcascade *cs, const AliESDEvent *esd) {
286
287 if (cs->Pt() < pMin) return kFALSE;
288
289 Int_t bidx=TMath::Abs(cs->GetBindex());
290 AliESDtrack *btrack=esd->GetTrack(bidx);
291 if (!AcceptTrack(btrack)) return kFALSE;
292
293 Float_t xy,z0;
294 btrack->GetImpactParameters(xy,z0);
295 if (TMath::Abs(xy)<0.03) return kFALSE;
296
297 const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
298 if (!vtx->GetStatus()) {
299 vtx=esd->GetPrimaryVertexTracks();
300 if (!vtx->GetStatus()) return kFALSE;
301 }
302 Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
303 if (cs->GetCascadeCosineOfPointingAngle(xv,yv,zv) < 0.999) return kFALSE;
304
305 if (cs->GetDcaXiDaughters() > 0.3) return kFALSE;
306
307 return kTRUE;
308}
309
1c102c9a 310static Bool_t AcceptPID(const AliPIDResponse *pidResponse,
311 const AliESDtrack *ptrack, AliStack *stack) {
312
f6357875 313 const AliExternalTrackParam *par=ptrack->GetInnerParam();
314 if (!par) return kTRUE;
315 if (par->GetP() > 1.) return kTRUE;
1c102c9a 316
317 if (stack) {
318 // MC PID
319 Int_t ntrk=stack->GetNtrack();
320 Int_t plab=TMath::Abs(ptrack->GetLabel());
321 if (plab>=0)
322 if (plab<ntrk) {
323 TParticle *pp=stack->Particle(plab);
c95561aa 324 if (pp->GetPDG()->Charge() > 0) {
f6357875 325 if (pp->GetPdgCode() == kProton) return kTRUE;
c95561aa 326 } else {
f6357875 327 if (pp->GetPdgCode() == kProtonBar) return kTRUE;
c95561aa 328 }
1c102c9a 329 }
330 } else {
331 // Real PID
f6357875 332 Double_t nsig=pidResponse->NumberOfSigmasTPC(ptrack,AliPID::kProton);
333 if (TMath::Abs(nsig) < 3.) return kTRUE;
1c102c9a 334 }
335
f6357875 336 return kFALSE;
1c102c9a 337}
338
e26966dc 339void AliAnalysisTaskCTauPbPb::UserExec(Option_t *)
340{
341
342 AliESDEvent *esd=(AliESDEvent *)InputEvent();
343
344 if (!esd) {
345 Printf("ERROR: esd not available");
346 return;
347 }
348
0b384e23 349 // Vertex selection
350 const AliESDVertex *vtx=esd->GetPrimaryVertexSPD();
351 if (!vtx->GetStatus()) {
352 vtx=esd->GetPrimaryVertexTracks();
353 if (!vtx->GetStatus()) return;
354 }
355 Double_t xv=vtx->GetXv(), yv=vtx->GetYv(), zv=vtx->GetZv();
356
357 if (TMath::Abs(zv) > 10.) return ;
358
e26966dc 359
360 // Physics selection
361 AliAnalysisManager *mgr= AliAnalysisManager::GetAnalysisManager();
362 AliInputEventHandler *hdr=(AliInputEventHandler*)mgr->GetInputEventHandler();
363 UInt_t maskIsSelected = hdr->IsEventSelected();
364 Bool_t isSelected = (maskIsSelected & AliVEvent::kMB);
365 if (!isSelected) return;
366
0b384e23 367
368 fMult->Fill(-100); //event counter
369
370
e26966dc 371 // Centrality selection
372 AliCentrality *cent=esd->GetCentrality();
373 if (!cent->IsEventInCentralityClass(fCMin,fCMax,"V0M")) return;
374
e26966dc 375 AliPIDResponse *pidResponse = hdr->GetPIDResponse();
376
e26966dc 377
378 //+++++++ MC
379 AliStack *stack = 0x0;
380 Double_t mcXv=0., mcYv=0., mcZv=0.;
381
382 if (fIsMC) {
383 AliMCEvent *mcEvent = MCEvent();
384 stack = mcEvent->Stack();
385 if (!stack) {
386 Printf("ERROR: stack not available");
387 return;
388 }
389
390 const AliVVertex *mcVtx=mcEvent->GetPrimaryVertex();
391
392 mcXv=mcVtx->GetX(); mcYv=mcVtx->GetY(); mcZv=mcVtx->GetZ();
393
394 Int_t ntrk=stack->GetNtrack(), ntrk0=ntrk;
395 while (ntrk--) {
396 TParticle *p0=stack->Particle(ntrk);
397 Int_t code=p0->GetPdgCode();
398 if (code != kK0Short)
0b384e23 399 if (code != kLambda0)
400 if (code != kLambda0Bar) continue;
e26966dc 401
402 Int_t plab=p0->GetFirstDaughter(), nlab=p0->GetLastDaughter();
403 if (nlab==plab) continue;
404 if (nlab<0) continue;
405 if (plab<0) continue;
406 if (nlab>=ntrk0) continue;
407 if (plab>=ntrk0) continue;
408 TParticle *part = stack->Particle(plab);
409 if (!part) continue;
410 TParticlePDG *partPDG = part->GetPDG();
411 if (!partPDG) continue;
412 Double_t charge=partPDG->Charge();
413 if (charge == 0.) continue;
414
415 Double_t pt=p0->Pt();
416 if (pt<pMin) continue;
417 if (TMath::Abs(p0->Y())>yMax) continue;
418
419 Double_t x=p0->Vx(), y=p0->Vy(), z=p0->Vz();
420 Double_t dx=mcXv-x, dy=mcYv-y, dz=mcZv-z;
421 Double_t l=TMath::Sqrt(dx*dx + dy*dy + dz*dz);
422
0b384e23 423 if (l > 0.001) continue; // secondary V0
e26966dc 424
425 x=part->Vx(); y=part->Vy();
426 dx=mcXv-x; dy=mcYv-y;
427 Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
428
429 if (code == kK0Short) {
430 fK0sMC->Fill(pt,lt);
431 }
432 if (code == kLambda0) {
433 fLambdaMC->Fill(pt,lt);
434 }
0b384e23 435 if (code == kLambda0Bar) {
436 fLambdaBarMC->Fill(pt,lt);
437 }
e26966dc 438 }
439 }
0b384e23 440 //+++++++
e26966dc 441
442
076bd7b1 443 Int_t ntrk1=esd->GetNumberOfTracks();
e26966dc 444 Int_t mult=0;
076bd7b1 445 for (Int_t i=0; i<ntrk1; i++) {
e26966dc 446 AliESDtrack *t=esd->GetTrack(i);
447 if (!t->IsOn(AliESDtrack::kTPCrefit)) continue;
448 Float_t xy,z0;
449 t->GetImpactParameters(xy,z0);
450 if (TMath::Abs(xy)>3.) continue;
451 if (TMath::Abs(z0)>3.) continue;
452 Double_t pt=t->Pt(),pz=t->Pz();
453 if (TMath::Abs(pz/pt)>0.8) continue;
454 mult++;
455
456 Double_t p=t->GetInnerParam()->GetP();
457 Double_t dedx=t->GetTPCsignal()/47.;
458 fdEdx->Fill(p,dedx,1);
459
1c102c9a 460 Double_t nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
e26966dc 461 if (TMath::Abs(nsig) < 3.) fdEdxPid->Fill(p,dedx,1);
462
463 }
464 fMult->Fill(mult);
465
466
467 Int_t nv0 = esd->GetNumberOfV0s();
468 while (nv0--) {
469 AliESDv0 *v0=esd->GetV0(nv0);
470
471 if (!AcceptV0(v0,esd)) continue;
472
473 Int_t nidx=TMath::Abs(v0->GetNindex());
474 AliESDtrack *ntrack=esd->GetTrack(nidx);
475 Int_t pidx=TMath::Abs(v0->GetPindex());
476 AliESDtrack *ptrack=esd->GetTrack(pidx);
477
478 Double_t x,y,z; v0->GetXYZ(x,y,z);
076bd7b1 479 Double_t dx1=x-xv, dy1=y-yv;
480 Double_t lt=TMath::Sqrt(dx1*dx1 + dy1*dy1);
e26966dc 481
482 Double_t pt=v0->Pt();
483
484 Bool_t ctK=kTRUE; if (0.4977*lt/pt > 3*2.68) ctK=kFALSE;
485 Bool_t ctL=kTRUE; if (1.1157*lt/pt > 3*7.89) ctL=kFALSE;
486
0b384e23 487 Bool_t isProton =AcceptPID(pidResponse, ptrack, stack);
488 Bool_t isProtonBar=AcceptPID(pidResponse, ntrack, stack);
1c102c9a 489
e26966dc 490 //+++++++ MC
491 if (stack) {
492 Int_t ntrk=stack->GetNtrack();
493
494 Int_t nlab=TMath::Abs(ntrack->GetLabel());
e26966dc 495 if (nlab<0) goto noas;
496 if (nlab>=ntrk) goto noas;
1c102c9a 497 TParticle *np=stack->Particle(nlab);
498
499 Int_t plab=TMath::Abs(ptrack->GetLabel());
e26966dc 500 if (plab<0) goto noas;
501 if (plab>=ntrk) goto noas;
e26966dc 502 TParticle *pp=stack->Particle(plab);
1c102c9a 503
e26966dc 504 Int_t i0=pp->GetFirstMother();
505 //if (np->GetFirstMother() != i0) goto noas;
506
507 Int_t in0=np->GetFirstMother();
508 if (in0<0) goto noas;
509 if (in0>=ntrk) goto noas;
510 if (in0 != i0) { // did the negative daughter decay ?
511 TParticle *nnp=stack->Particle(in0);
512 if (nnp->GetFirstMother() != i0) goto noas;
513 }
514
515 if (i0<0) goto noas;
516 if (i0>=ntrk) goto noas;
517 TParticle *p0=stack->Particle(i0);
518
519 Int_t code=p0->GetPdgCode();
520 if (code != kK0Short)
0b384e23 521 if (code != kLambda0)
522 if (code != kLambda0Bar) goto noas;
e26966dc 523
524 if (p0->Pt()<pMin) goto noas;
525 if (TMath::Abs(p0->Y())>yMax ) goto noas;
526
527
528 Double_t dz=mcZv - p0->Vz(), dy=mcYv - p0->Vy(), dx=mcXv - p0->Vx();
529 Double_t l = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
530
531 dx = mcXv - pp->Vx(); dy = mcYv - pp->Vy();
532 Double_t ltAs=TMath::Sqrt(dx*dx + dy*dy);
533 Double_t ptAs=p0->Pt();
534
0b384e23 535 if (l > 0.001) { // Secondary V0
536 if (code != kLambda0)
537 if (code != kLambda0Bar) goto noas;
e26966dc 538 Int_t nx=p0->GetFirstMother();
539 if (nx<0) goto noas;
540 if (nx>=ntrk) goto noas;
541 TParticle *xi=stack->Particle(nx);
542 Int_t xcode=xi->GetPdgCode();
0b384e23 543 if (code == kLambda0) {
544 if ( xcode != kXiMinus )
545 if ( xcode != 3322 ) goto noas;
f6357875 546 if (ctL) fLambdaFromXi->Fill(ptAs,ltAs,xi->Pt());
0b384e23 547 } else if (code == kLambda0Bar) {
548 if ( xcode != kXiPlusBar )
549 if ( xcode != -3322 ) goto noas;
f6357875 550 if (ctL) fLambdaBarFromXiBar->Fill(ptAs,ltAs,xi->Pt());
0b384e23 551 }
e26966dc 552 } else {
553 if (code == kLambda0) {
554 if (ctL) fLambdaAs->Fill(ptAs,ltAs);
0b384e23 555 } else if (code == kLambda0Bar) {
556 if (ctL) fLambdaBarAs->Fill(ptAs,ltAs);
e26966dc 557 } else {
558 if (ctK) fK0sAs->Fill(ptAs,ltAs);
559 }
560 }
561
562 }
563 //++++++++
564
565 noas:
566
567 Double_t dca=v0->GetDcaV0Daughters();
568 Double_t cpa=v0->GetV0CosineOfPointingAngle();
569
570 Double_t mass=0., m=0., s=0.;
571 if (ctK)
572 if (TMath::Abs(v0->RapK0Short())<yMax) {
573 v0->ChangeMassHypothesis(kK0Short);
574
575 mass=v0->GetEffMass();
576 fK0sM->Fill(mass,pt);
577
578 m=TDatabasePDG::Instance()->GetParticle(kK0Short)->Mass();
579 s=0.0044 + (0.008-0.0044)/(10-1)*(pt - 1.);
580 if (TMath::Abs(m-mass) < 3*s) {
581 fK0sSi->Fill(pt,lt);
582 }
583 if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
584 fK0sSi->Fill(pt,lt,-1);
585 }
586 if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
587 fK0sSi->Fill(pt,lt,-1);
588 }
589 }
590
591 if (ctL)
1c102c9a 592 if (isProton)
e26966dc 593 if (TMath::Abs(v0->RapLambda())<yMax) {
e26966dc 594 v0->ChangeMassHypothesis(kLambda0);
595
596 mass=v0->GetEffMass();
597 fLambdaM->Fill(mass,pt);
598
599 m=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();
600 //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
601 //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
602 s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
603 if (TMath::Abs(m-mass) < 3*s) {
604 fLambdaSi->Fill(pt,lt);
605 fCPA->Fill(cpa,1);
606 fDCA->Fill(dca,1);
607 }
608 if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
609 fLambdaSi->Fill(pt,lt,-1);
610 fCPA->Fill(cpa,-1);
611 fDCA->Fill(dca,-1);
612 }
613 if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
614 fLambdaSi->Fill(pt,lt,-1);
615 fCPA->Fill(cpa,-1);
616 fDCA->Fill(dca,-1);
617 }
618 }
0b384e23 619
620 if (ctL)
621 if (isProtonBar)
622 if (TMath::Abs(v0->RapLambda())<yMax) {
623 v0->ChangeMassHypothesis(kLambda0Bar);
624
625 mass=v0->GetEffMass();
626 fLambdaBarM->Fill(mass,pt);
627
628 m=TDatabasePDG::Instance()->GetParticle(kLambda0Bar)->Mass();
629 //s=0.0027 + (0.004-0.0027)/(10-1)*(pt-1);
630 //s=0.0015 + (0.002-0.0015)/(2.6-1)*(pt-1);
631 s=0.0023 + (0.004-0.0023)/(6-1)*(pt-1);
632 if (TMath::Abs(m-mass) < 3*s) {
633 fLambdaBarSi->Fill(pt,lt);
634 fCPA->Fill(cpa,1);
635 fDCA->Fill(dca,1);
636 }
637 if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
638 fLambdaBarSi->Fill(pt,lt,-1);
639 fCPA->Fill(cpa,-1);
640 fDCA->Fill(dca,-1);
641 }
642 if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
643 fLambdaBarSi->Fill(pt,lt,-1);
644 fCPA->Fill(cpa,-1);
645 fDCA->Fill(dca,-1);
646 }
647 }
e26966dc 648 }
649
650 Double_t kine0;
651 Int_t ncs=esd->GetNumberOfCascades();
652 for (Int_t i=0; i<ncs; i++) {
653 AliESDcascade *cs=esd->GetCascade(i);
654
655 if (TMath::Abs(cs->RapXi()) > yMax) continue;
656 if (!AcceptCascade(cs,esd)) continue;
657
658 AliESDv0 *v0 = (AliESDv0*)cs;
659 if (TMath::Abs(v0->RapLambda()) > yMax) continue;
660 if (!AcceptV0(v0,esd)) continue;
661
662 Double_t pt=cs->Pt();
663
1c102c9a 664 Int_t pidx=TMath::Abs(v0->GetPindex());
665 AliESDtrack *ptrack=esd->GetTrack(pidx);
0b384e23 666 Bool_t isProton =AcceptPID(pidResponse, ptrack, stack);
1c102c9a 667
0b384e23 668 Int_t nidx=TMath::Abs(v0->GetNindex());
669 AliESDtrack *ntrack=esd->GetTrack(nidx);
670 Bool_t isProtonBar=AcceptPID(pidResponse, ntrack, stack);
1c102c9a 671
e26966dc 672 Int_t charge=cs->Charge();
1c102c9a 673 if (isProton)
e26966dc 674 if (charge < 0) {
e26966dc 675 cs->ChangeMassHypothesis(kine0,kXiMinus);
676 Double_t mass=cs->GetEffMassXi();
677 pt=cs->Pt();
678 fXiM->Fill(mass,pt);
679 Double_t m=TDatabasePDG::Instance()->GetParticle(kXiMinus)->Mass();
680 //Double_t s=0.0037;
681 Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
682 if (TMath::Abs(m-mass) < 3*s) {
683 fXiSiP->Fill(pt);
684 }
685 if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
686 fXiSiP->Fill(pt,-1);
687 }
688 if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
689 fXiSiP->Fill(pt,-1);
690 }
691 }
0b384e23 692 if (isProtonBar)
693 if (charge > 0) {
694 cs->ChangeMassHypothesis(kine0,kXiPlusBar);
695 Double_t mass=cs->GetEffMassXi();
696 pt=cs->Pt();
697 fXiBarM->Fill(mass,pt);
698 Double_t m=TDatabasePDG::Instance()->GetParticle(kXiPlusBar)->Mass();
699 //Double_t s=0.0037;
700 Double_t s=0.002 + (0.0032-0.002)/(6-1.5)*(pt-1.5);
701 if (TMath::Abs(m-mass) < 3*s) {
702 fXiBarSiP->Fill(pt);
703 }
704 if (TMath::Abs(m-mass + 4.5*s) < 1.5*s) {
705 fXiBarSiP->Fill(pt,-1);
706 }
707 if (TMath::Abs(m-mass - 4.5*s) < 1.5*s) {
708 fXiBarSiP->Fill(pt,-1);
709 }
710 }
e26966dc 711 }
712
713}
714
715void AliAnalysisTaskCTauPbPb::Terminate(Option_t *)
716{
717 // The Terminate() function is the last function to be called during
718 // a query. It always runs on the client, it can be used to present
719 // the results graphically or save the results to file.
720
721 fOutput=(TList*)GetOutputData(1);
722 if (!fOutput) {
723 Printf("ERROR: fOutput not available");
724 return;
725 }
726
727 fMult = dynamic_cast<TH1F*>(fOutput->FindObject("fMult")) ;
728 if (!fMult) {
729 Printf("ERROR: fMult not available");
730 return;
731 }
732
733 fdEdx = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdx")) ;
734 if (!fdEdx) {
735 Printf("ERROR: fdEdx not available");
736 return;
737 }
738
739 fdEdxPid = dynamic_cast<TH2F*>(fOutput->FindObject("fdEdxPid")) ;
740 if (!fdEdxPid) {
741 Printf("ERROR: fdEdxPid not available");
742 return;
743 }
744
745
746 fK0sMC = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sMC")) ;
747 if (!fK0sMC) {
748 Printf("ERROR: fK0sMC not available");
749 return;
750 }
751 TH1D *k0sMcPx=fK0sMC->ProjectionX(); k0sMcPx->Sumw2();
752 fK0sAs = dynamic_cast<TH2F*>(fOutput->FindObject("fK0sAs")) ;
753 if (!fK0sAs) {
754 Printf("ERROR: fK0sAs not available");
755 return;
756 }
757 TH1D *k0sAsPx=fK0sAs->ProjectionX();
758 k0sAsPx->Sumw2(); //k0sAsPx->Scale(0.69);
759
760
761
762 fLambdaFromXi = dynamic_cast<TH3F*>(fOutput->FindObject("fLambdaFromXi")) ;
763 if (!fLambdaFromXi) {
764 Printf("ERROR: fLambdaFromXi not available");
765 return;
766 }
767 TH1D *lambdaFromXiPx=fLambdaFromXi->ProjectionX(); lambdaFromXiPx->Sumw2();
768
769
770 fLambdaMC = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaMC")) ;
771 if (!fLambdaMC) {
772 Printf("ERROR: fLambdaMC not available");
773 return;
774 }
775 TH1D *lambdaMcPx=fLambdaMC->ProjectionX(); lambdaMcPx->Sumw2();
776
777 fLambdaAs = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaAs")) ;
778 if (!fLambdaAs) {
779 Printf("ERROR: fLambdaAs not available");
780 return;
781 }
782 TH1D *lambdaAsPx=fLambdaAs->ProjectionX();
783 lambdaAsPx->Sumw2(); //lambdaAsPx->Scale(0.64);
784
785 fLambdaSi = dynamic_cast<TH2F*>(fOutput->FindObject("fLambdaSi")) ;
786 if (!fLambdaSi) {
787 Printf("ERROR: fLambdaSi not available");
788 return;
789 }
790 TH1D *lambdaSiPx=fLambdaSi->ProjectionX();
791 lambdaSiPx->SetName("fLambdaPt");
792 lambdaSiPx->Sumw2();
793
794 fLambdaEff = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaEff")) ;
795 if (!fLambdaEff) {
796 Printf("ERROR: fLambdaEff not available");
797 return;
798 }
799 fLambdaPt = dynamic_cast<TH1D*>(fOutput->FindObject("fLambdaPt")) ;
800 if (!fLambdaPt) {
801 Printf("ERROR: fLambdaPt not available");
802 return;
803 }
804
805
806 if (!gROOT->IsBatch()) {
807
808 TCanvas *c1 = new TCanvas("c1","Mulitplicity");
809 c1->SetLogy();
810 fMult->DrawCopy() ;
811
812 new TCanvas("c2","dE/dx");
813 fdEdx->DrawCopy() ;
814
815 new TCanvas("c3","dE/dx with PID");
816 fdEdxPid->DrawCopy() ;
817
818 if (fIsMC) {
819 /*
820 TH1D effK(*k0sAsPx); effK.SetTitle("Efficiency for K0s");
821 effK.Divide(k0sAsPx,k0sMcPx,1,1,"b");
822 new TCanvas("c4","Efficiency for K0s");
823 effK.DrawCopy("E") ;
824 */
825
826 fLambdaEff->Divide(lambdaAsPx,lambdaMcPx,1,1,"b");
827 new TCanvas("c5","Efficiency for #Lambda");
828 fLambdaEff->DrawCopy("E") ;
829
830 lambdaSiPx->Add(lambdaFromXiPx,-1);
831 lambdaSiPx->Divide(fLambdaEff);
832
833 new TCanvas("c6","Corrected #Lambda pt");
834 lambdaSiPx->SetTitle("Corrected #Lambda pt");
835 *fLambdaPt = *lambdaSiPx;
836 fLambdaPt->SetLineColor(2);
837 fLambdaPt->DrawCopy("E");
838
839 lambdaMcPx->DrawCopy("same");
840
841 } else {
842 new TCanvas("c6","Raw #Lambda pt");
843 lambdaSiPx->SetTitle("Raw #Lambda pt");
844 *fLambdaPt = *lambdaSiPx;
845 fLambdaPt->SetLineColor(2);
846 fLambdaPt->DrawCopy("E");
847 }
848 }
849}