]>
Commit | Line | Data |
---|---|---|
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 | ||
30 | extern TROOT *gROOT; | |
31 | ||
32 | ClassImp(AliAnalysisTaskCTauPbPb) | |
33 | ||
25aeb948 | 34 | static Int_t nbins=100; // number of bins |
e26966dc | 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 | ||
0b384e23 | 64 | fLambdaBarM(0), |
65 | fLambdaBarSi(0), | |
66 | fLambdaBarMC(0), | |
67 | fLambdaBarAs(0), | |
68 | ||
e26966dc | 69 | fCPA(0), |
70 | fDCA(0), | |
71 | ||
72 | fLambdaEff(0), | |
73 | fLambdaPt(0), | |
74 | ||
75 | fLambdaFromXi(0), | |
76 | fXiM(0), | |
0b384e23 | 77 | fXiSiP(0), |
78 | ||
79 | fLambdaBarFromXiBar(0), | |
80 | fXiBarM(0), | |
81 | fXiBarSiP(0) | |
e26966dc | 82 | { |
83 | // Constructor. Initialization of pointers | |
84 | DefineOutput(1, TList::Class()); | |
85 | } | |
86 | ||
87 | void 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 | ||
232 | static 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 | ||
247 | static 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 | ||
285 | static 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 | 310 | static 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 | 339 | void 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 | ||
715 | void 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 | } |