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