Using the "perfect PID" when running on Monte Carlo
[u/mrichter/AliRoot.git] / PWGLF / STRANGENESS / LambdaK0PbPb / AliAnalysisTaskCTauPbPb.cxx
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
30 extern TROOT *gROOT;
31
32 ClassImp(AliAnalysisTaskCTauPbPb)
33
34 static Int_t    nbins=102;  // number of bins
35 static Double_t lMin=0.0, lMax=100.;
36 static Double_t pMin=0.0, pMax=10.;
37 static 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
44 AliAnalysisTaskCTauPbPb::AliAnalysisTaskCTauPbPb(const char *name) :
45 AliAnalysisTaskSE(name),
46 fIsMC(kFALSE),
47 fCMin(0.),
48 fCMax(90.),
49 fOutput(0),
50 fMult(0),
51 fdEdx(0),
52 fdEdxPid(0),
53
54 fK0sM(0),
55 fK0sSi(0),
56 fK0sMC(0),
57 fK0sAs(0),
58
59 fLambdaM(0),
60 fLambdaSi(0),
61 fLambdaMC(0),
62 fLambdaAs(0),
63
64 fCPA(0),
65 fDCA(0),
66
67 fLambdaEff(0),
68 fLambdaPt(0),
69
70 fLambdaFromXi(0),
71 fXiM(0),
72 fXiSiP(0)
73 {
74   // Constructor. Initialization of pointers
75   DefineOutput(1, TList::Class());
76 }
77
78 void 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
182 static 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
192   if (TMath::Abs(t->Eta()) > 0.8) return kFALSE;
193
194   return kTRUE;   
195 }
196
197 static 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
235 static 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
260 static 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
287 void 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
383   Int_t ntrk=esd->GetNumberOfTracks();
384   Int_t mult=0;
385   for (Int_t i=0; i<ntrk; i++) {
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
400     Double_t nsig=pidResponse->NumberOfSigmasTPC(t,AliPID::kProton);
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);
419       Double_t dx=x-xv, dy=y-yv;
420       Double_t lt=TMath::Sqrt(dx*dx + dy*dy);
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
427       Bool_t isProton=AcceptPID(pidResponse, ptrack, stack);
428
429       //+++++++ MC
430       if (stack) {
431          Int_t ntrk=stack->GetNtrack();
432
433          Int_t nlab=TMath::Abs(ntrack->GetLabel());
434          if (nlab<0) goto noas;      
435          if (nlab>=ntrk) goto noas;      
436          TParticle *np=stack->Particle(nlab);
437
438          Int_t plab=TMath::Abs(ptrack->GetLabel());
439          if (plab<0) goto noas;      
440          if (plab>=ntrk) goto noas;      
441          TParticle *pp=stack->Particle(plab);
442
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)
521       if (isProton)
522       if (TMath::Abs(v0->RapLambda())<yMax) {
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
564       Int_t pidx=TMath::Abs(v0->GetPindex());
565       AliESDtrack *ptrack=esd->GetTrack(pidx);
566
567       Bool_t isProton=AcceptPID(pidResponse, ptrack, stack);
568
569       Int_t charge=cs->Charge();      
570       if (isProton)
571       if (charge < 0) {         
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
593 void 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 }