]>
Commit | Line | Data |
---|---|---|
10757ee9 | 1 | |
2 | #include <TROOT.h> | |
3 | #include <TChain.h> | |
4 | #include <TFile.h> | |
5 | #include <TH3F.h> | |
6 | // | |
7 | #include "TParticle.h" | |
8 | #include "TDatabasePDG.h" | |
9 | #include "AliRunLoader.h" | |
10 | #include "AliStack.h" | |
11 | ||
12 | ||
13 | ||
14 | #include <TPDGCode.h> | |
15 | #include <TStyle.h> | |
16 | #include "TLinearFitter.h" | |
17 | #include "TMatrixD.h" | |
18 | #include "TTreeStream.h" | |
19 | #include "TF1.h" | |
20 | ||
21 | ||
22 | ||
23 | #include "AliMagF.h" | |
24 | #include "AliTracker.h" | |
25 | #include "AliESD.h" | |
26 | #include "AliESDtrack.h" | |
27 | #include "AliESDfriend.h" | |
28 | #include "AliESDfriendTrack.h" | |
29 | #include "AliTPCseed.h" | |
30 | #include "AliTPCclusterMI.h" | |
31 | ||
32 | #include "AliKFParticle.h" | |
33 | #include "AliKFVertex.h" | |
34 | ||
35 | #include "AliTrackPointArray.h" | |
36 | #include "TCint.h" | |
37 | #include "AliTPCcalibV0.h" | |
38 | #include "AliV0.h" | |
39 | ||
40 | ||
41 | ||
42 | ||
43 | ||
44 | ClassImp(AliTPCcalibV0) | |
45 | ||
46 | ||
47 | AliTPCcalibV0::AliTPCcalibV0() : | |
48 | TNamed(), | |
49 | fDebugStream(0), | |
50 | fStack(0), | |
51 | fOutput(0), | |
52 | fESD(0), | |
53 | fPdg(0), | |
54 | fParticles(0), | |
55 | fV0s(0), | |
56 | fGammas(0), | |
57 | fV0type(0), | |
58 | fTPCdEdx(0), | |
59 | fTPCdEdxPi(0), | |
60 | fTPCdEdxEl(0), | |
61 | fTPCdEdxP(0) | |
62 | { | |
63 | G__SetCatchException(0); | |
64 | fDebugStream = new TTreeSRedirector("V0debug.root"); | |
65 | fPdg = new TDatabasePDG; | |
66 | ||
67 | ||
68 | // create output histograms | |
69 | fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); | |
70 | BinLogX(fTPCdEdx); | |
71 | fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); | |
72 | fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); | |
73 | fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300); | |
74 | ||
75 | ||
76 | } | |
77 | ||
78 | AliTPCcalibV0::~AliTPCcalibV0(){ | |
79 | // | |
80 | // | |
81 | // | |
82 | delete fDebugStream; | |
83 | } | |
84 | ||
85 | ||
86 | ||
87 | void AliTPCcalibV0::ProofSlaveBegin(TList * output) | |
88 | { | |
89 | // Called on PROOF - fill output list | |
90 | } | |
91 | ||
92 | ||
93 | void AliTPCcalibV0::ProcessESD(AliESD *esd, AliStack *stack){ | |
94 | // | |
95 | // | |
96 | // | |
97 | fESD = esd; | |
98 | AliKFParticle::SetField(esd->GetMagneticField()); | |
99 | MakeV0s(); | |
100 | if (stack) { | |
101 | fStack = stack; | |
102 | MakeMC(); | |
103 | }else{ | |
104 | fStack =0; | |
105 | } | |
106 | } | |
107 | ||
108 | void AliTPCcalibV0::MakeMC(){ | |
109 | // | |
110 | // MC comparison | |
111 | // 1. Select interesting particles | |
112 | // 2. Assign the recosntructed particles | |
113 | // | |
114 | //1. Select interesting particles | |
115 | const Float_t kMinP = 0.2; | |
116 | const Float_t kMinPt = 0.1; | |
117 | const Float_t kMaxR = 0.5; | |
118 | const Float_t kMaxTan = 1.2; | |
119 | const Float_t kMaxRad = 150; | |
120 | // | |
121 | if (!fParticles) fParticles = new TObjArray; | |
122 | TParticle *part=0; | |
123 | // | |
124 | Int_t entries = fStack->GetNtrack(); | |
125 | for (Int_t ipart=0; ipart<entries; ipart++){ | |
126 | part = fStack->Particle(ipart); | |
127 | if (!part) continue; | |
128 | if (part->P()<kMinP) continue; | |
129 | if (part->R()>kMaxR) continue; | |
130 | if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue; | |
131 | Bool_t isInteresting = kFALSE; | |
132 | if (part->GetPdgCode()==22) isInteresting =kTRUE; | |
133 | if (part->GetPdgCode()==310) isInteresting =kTRUE; | |
134 | if (part->GetPdgCode()==111) isInteresting =kTRUE; | |
135 | if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE; | |
136 | ||
137 | // | |
138 | if (!isInteresting) continue; | |
139 | fParticles->AddLast(new TParticle(*part)); | |
140 | } | |
141 | if (fParticles->GetEntries()<1) { | |
142 | return; | |
143 | } | |
144 | // | |
145 | // | |
146 | // | |
147 | Int_t sentries=fParticles->GetEntries();; | |
148 | for (Int_t ipart=0; ipart<sentries; ipart++){ | |
149 | TParticle *part = (TParticle*)fParticles->At(ipart); | |
150 | TParticle *p0 = 0; | |
151 | TParticle *p1 = 0; | |
152 | ||
153 | Int_t nold =0; | |
154 | Int_t nnew =0; | |
155 | Int_t id0 = part->GetDaughter(0); | |
156 | Int_t id1 = part->GetDaughter(1); | |
157 | if (id0>=fStack->GetNtrack() ) id0*=-1; | |
158 | if (id1>=fStack->GetNtrack() ) id1*=-1; | |
159 | Bool_t findable = kTRUE; | |
160 | if (id0<0 || id1<0) findable = kFALSE; | |
161 | Int_t charge =0; | |
162 | if (findable){ | |
163 | p0 = fStack->Particle(id0); | |
164 | if (p0->R()>kMaxRad) findable = kFALSE; | |
165 | if (p0->Pt()<kMinPt) findable = kFALSE; | |
166 | if (p0->Vz()>250) findable= kFALSE; | |
167 | if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; | |
168 | if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE; | |
169 | else | |
170 | if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++; | |
171 | ||
172 | p1 = fStack->Particle(id1); | |
173 | if (p1->R()>kMaxRad) findable = kFALSE; | |
174 | if (p1->Pt()<kMinPt) findable = kFALSE; | |
175 | if (TMath::Abs(p1->Vz())>250) findable= kFALSE; | |
176 | if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE; | |
177 | if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE; | |
178 | else | |
179 | if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++; | |
180 | ||
181 | } | |
182 | // (*fDebugStream)<<"MC0"<< | |
183 | // "P.="<<part<< | |
184 | // "findable="<<findable<< | |
185 | // "id0="<<id0<< | |
186 | // "id1="<<id1<< | |
187 | // "\n"; | |
188 | if (!findable) continue; | |
189 | Float_t minpt = TMath::Min(p0->Pt(), p1->Pt()); | |
190 | Int_t type=-1; | |
191 | ||
192 | // | |
193 | // | |
194 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
195 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
196 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
197 | // select coresponding track | |
198 | AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); | |
199 | if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue; | |
200 | AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); | |
201 | if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue; | |
202 | TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel())); | |
203 | TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel())); | |
204 | // | |
205 | // | |
206 | if ( v0->GetOnFlyStatus()) nnew++; | |
207 | if (!v0->GetOnFlyStatus()) nold++; | |
208 | if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 && TMath::Abs(pp->GetPdgCode())==11) | |
209 | type =1; | |
210 | if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 && TMath::Abs(pp->GetPdgCode())==211) | |
211 | type =0; | |
212 | if (part->GetPdgCode()==3122){ | |
213 | if (TMath::Abs(pn->GetPdgCode())==210 ) type=2; | |
214 | else type=3; | |
215 | } | |
216 | AliKFParticle *v0kf = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); | |
217 | v0kf->SetProductionVertex( primVtx ); | |
218 | AliKFParticle *v0kfc = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode()); | |
219 | v0kfc->SetProductionVertex( primVtx ); | |
220 | v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass()); | |
221 | Float_t chi2 = v0kf->GetChi2(); | |
222 | Float_t chi2C = v0kf->GetChi2(); | |
223 | // | |
224 | // | |
225 | (*fDebugStream)<<"MCRC"<< | |
226 | "P.="<<part<< | |
227 | "type="<<type<< | |
228 | "chi2="<<chi2<< | |
229 | "chi2C="<<chi2C<< | |
230 | "minpt="<<minpt<< | |
231 | "id0="<<id0<< | |
232 | "id1="<<id1<< | |
233 | "Pn.="<<pn<< | |
234 | "Pp.="<<pp<< | |
235 | "tn.="<<trackN<< | |
236 | "tp.="<<trackP<< | |
237 | "nold.="<<nold<< | |
238 | "nnew.="<<nnew<< | |
239 | "v0.="<<v0<< | |
240 | "v0kf.="<<v0kf<< | |
241 | "v0kfc.="<<v0kfc<< | |
242 | "\n"; | |
243 | delete v0kf; | |
244 | delete v0kfc; | |
245 | // | |
246 | } | |
247 | (*fDebugStream)<<"MC"<< | |
248 | "P.="<<part<< | |
249 | "charge="<<charge<< | |
250 | "type="<<type<< | |
251 | "minpt="<<minpt<< | |
252 | "id0="<<id0<< | |
253 | "id1="<<id1<< | |
254 | "P0.="<<p0<< | |
255 | "P1.="<<p1<< | |
256 | "nold="<<nold<< | |
257 | "nnew="<<nnew<< | |
258 | "\n"; | |
259 | } | |
260 | fParticles->Delete(); | |
261 | ||
262 | } | |
263 | ||
264 | ||
265 | ||
266 | void AliTPCcalibV0::MakeV0s(){ | |
267 | // | |
268 | // | |
269 | // | |
270 | const Int_t kMinCluster=40; | |
271 | const Float_t kMinR =0; | |
272 | if (! fV0s) fV0s = new TObjArray(10); | |
273 | fV0s->Clear(); | |
274 | // | |
275 | // Old V0 finder | |
276 | // | |
277 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
278 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
279 | if (v0->GetOnFlyStatus()) continue; | |
280 | fV0s->AddLast(v0); | |
281 | } | |
282 | ProcessV0(0); | |
283 | fV0s->Clear(0); | |
284 | // | |
285 | // MI V0 finder | |
286 | // | |
287 | for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){ | |
288 | AliESDv0 * v0 = fESD->GetV0(ivertex); | |
289 | if (!v0->GetOnFlyStatus()) continue; | |
290 | fV0s->AddLast(v0); | |
291 | } | |
292 | ProcessV0(1); | |
293 | fV0s->Clear(); | |
294 | return; | |
295 | // | |
296 | // combinatorial | |
297 | // | |
298 | Int_t ntracks = fESD->GetNumberOfTracks(); | |
299 | for (Int_t itrack0=0; itrack0<ntracks; itrack0++){ | |
300 | AliESDtrack * track0 = fESD->GetTrack(itrack0); | |
301 | if (track0->GetSign()>0) continue; | |
302 | if ( track0->GetTPCNcls()<kMinCluster) continue; | |
303 | if (track0->GetKinkIndex(0)>0) continue; | |
304 | // | |
305 | for (Int_t itrack1=0; itrack1<ntracks; itrack1++){ | |
306 | AliESDtrack * track1 = fESD->GetTrack(itrack1); | |
307 | if (track1->GetSign()<0) continue; | |
308 | if ( track1->GetTPCNcls()<kMinCluster) continue; | |
309 | if (track1->GetKinkIndex(0)>0) continue; | |
310 | // | |
311 | // AliExternalTrackParam param0(*track0); | |
312 | // AliExternalTrackParam param1(*track1); | |
313 | AliV0 vertex; | |
314 | vertex.SetParamN(*track0); | |
315 | vertex.SetParamP(*track1); | |
316 | Float_t xyz[3]; | |
317 | xyz[0] = fESD->GetPrimaryVertex()->GetXv(); | |
318 | xyz[1] = fESD->GetPrimaryVertex()->GetYv(); | |
319 | xyz[2] = fESD->GetPrimaryVertex()->GetZv(); | |
320 | vertex.Update(xyz); | |
321 | if (vertex.GetRr()<kMinR) continue; | |
322 | if (vertex.GetDcaV0Daughters()>1.) continue; | |
323 | if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue; | |
324 | // if (vertex.GetPointAngle()<0.9) continue; | |
325 | vertex.SetIndex(0,itrack0); | |
326 | vertex.SetIndex(1,itrack1); | |
327 | fV0s->AddLast(new AliV0(vertex)); | |
328 | } | |
329 | } | |
330 | ProcessV0(2); | |
331 | for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i); | |
332 | fV0s->Clear(); | |
333 | } | |
334 | ||
335 | ||
336 | ||
337 | ||
338 | ||
339 | ||
340 | // void AliTPCcalibV0::ProcessV0(Int_t ftype){ | |
341 | // // | |
342 | // // | |
343 | // const Double_t ktimeK0 = 2.684; | |
344 | // const Double_t ktimeLambda = 7.89; | |
345 | ||
346 | ||
347 | // if (! fGammas) fGammas = new TObjArray(10); | |
348 | // fGammas->Clear(); | |
349 | // Int_t nV0s = fV0s->GetEntries(); | |
350 | // if (nV0s==0) return; | |
351 | // AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
352 | // // | |
353 | // for (Int_t ivertex=0; ivertex<nV0s; ivertex++){ | |
354 | // AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex); | |
355 | // AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); | |
356 | // AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); | |
357 | // // | |
358 | // // | |
359 | // // | |
360 | // AliKFParticle *v0K0 = Fit(primVtx,v0,211,211); | |
361 | // AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); | |
362 | // AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211); | |
363 | // AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212); | |
364 | // //Set production vertex | |
365 | // v0K0->SetProductionVertex( primVtx ); | |
366 | // v0Gamma->SetProductionVertex( primVtx ); | |
367 | // v0Lambda42->SetProductionVertex( primVtx ); | |
368 | // v0Lambda24->SetProductionVertex( primVtx ); | |
369 | // Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; | |
370 | // v0K0->GetMass( massK0,massSigma); | |
371 | // v0Gamma->GetMass( massGamma,massSigma); | |
372 | // v0Lambda42->GetMass( massLambda42,massSigma); | |
373 | // v0Lambda24->GetMass( massLambda24,massSigma); | |
374 | // Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); | |
375 | // Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); | |
376 | // Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); | |
377 | // Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); | |
378 | // // | |
379 | // // Mass Contrained params | |
380 | // // | |
381 | // AliKFParticle *v0K0C = Fit(primVtx,v0,211,211); | |
382 | // AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); | |
383 | // AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211); | |
384 | // AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212); | |
385 | // // | |
386 | // v0K0C->SetProductionVertex( primVtx ); | |
387 | // v0GammaC->SetProductionVertex( primVtx ); | |
388 | // v0Lambda42C->SetProductionVertex( primVtx ); | |
389 | // v0Lambda24C->SetProductionVertex( primVtx ); | |
390 | ||
391 | // v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); | |
392 | // v0GammaC->SetMassConstraint(0); | |
393 | // v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); | |
394 | // v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); | |
395 | // // | |
396 | // Double_t timeK0, sigmaTimeK0; | |
397 | // Double_t timeLambda42, sigmaTimeLambda42; | |
398 | // Double_t timeLambda24, sigmaTimeLambda24; | |
399 | // v0K0C->GetLifeTime(timeK0, sigmaTimeK0); | |
400 | // //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); | |
401 | // v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); | |
402 | // v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); | |
403 | ||
404 | ||
405 | // // | |
406 | // Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); | |
407 | // if (chi2K0C<0) chi2K0C=100; | |
408 | // Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); | |
409 | // if (chi2GammaC<0) chi2GammaC=100; | |
410 | // Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); | |
411 | // if (chi2Lambda42C<0) chi2Lambda42C=100; | |
412 | // Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); | |
413 | // if (chi2Lambda24C<0) chi2Lambda24C=100; | |
414 | // // | |
415 | // Float_t minChi2C=99; | |
416 | // Int_t type =-1; | |
417 | // if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;} | |
418 | // if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;} | |
419 | // if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;} | |
420 | // if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;} | |
421 | // Float_t minChi2=99; | |
422 | // Int_t type0 =-1; | |
423 | // if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;} | |
424 | // if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;} | |
425 | // if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;} | |
426 | // if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;} | |
427 | // Float_t betaGammaP = trackN->GetP()/fPdg->GetParticle(-2212)->Mass(); | |
428 | // Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass(); | |
429 | // Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass(); | |
430 | // Float_t dedxTeorP = TPCBetheBloch(betaGammaP); | |
431 | // Float_t dedxTeorPi = TPCBetheBloch(betaGammaPi);; | |
432 | // Float_t dedxTeorEl = TPCBetheBloch(betaGammaEl);; | |
433 | // // | |
434 | // // | |
435 | // if (minChi2>50) continue; | |
436 | // (*fDebugStream)<<"V0"<< | |
437 | // "ftype="<<ftype<< | |
438 | // "v0.="<<v0<< | |
439 | // "trackN.="<<trackN<< | |
440 | // "trackP.="<<trackP<< | |
441 | // // | |
442 | // "dedxTeorP="<<dedxTeorP<< | |
443 | // "dedxTeorPi="<<dedxTeorPi<< | |
444 | // "dedxTeorEl="<<dedxTeorEl<< | |
445 | // // | |
446 | // "type="<<type<< | |
447 | // "chi2C="<<minChi2C<< | |
448 | // "v0K0.="<<v0K0<< | |
449 | // "v0Gamma.="<<v0Gamma<< | |
450 | // "v0Lambda42.="<<v0Lambda42<< | |
451 | // "v0Lambda24.="<<v0Lambda24<< | |
452 | // // | |
453 | // "chi20K0.="<<chi2K0<< | |
454 | // "chi2Gamma.="<<chi2Gamma<< | |
455 | // "chi2Lambda42.="<<chi2Lambda42<< | |
456 | // "chi2Lambda24.="<<chi2Lambda24<< | |
457 | // // | |
458 | // "chi20K0c.="<<chi2K0C<< | |
459 | // "chi2Gammac.="<<chi2GammaC<< | |
460 | // "chi2Lambda42c.="<<chi2Lambda42C<< | |
461 | // "chi2Lambda24c.="<<chi2Lambda24C<< | |
462 | // // | |
463 | // "v0K0C.="<<v0K0C<< | |
464 | // "v0GammaC.="<<v0GammaC<< | |
465 | // "v0Lambda42C.="<<v0Lambda42C<< | |
466 | // "v0Lambda24C.="<<v0Lambda24C<< | |
467 | // // | |
468 | // "massK0="<<massK0<< | |
469 | // "massGamma="<<massGamma<< | |
470 | // "massLambda42="<<massLambda42<< | |
471 | // "massLambda24="<<massLambda24<< | |
472 | // // | |
473 | // "timeK0="<<timeK0<< | |
474 | // "timeLambda42="<<timeLambda42<< | |
475 | // "timeLambda24="<<timeLambda24<< | |
476 | // "\n"; | |
477 | // if (type==1) fGammas->AddLast(v0); | |
478 | // // | |
479 | // // | |
480 | // // | |
481 | // delete v0K0; | |
482 | // delete v0Gamma; | |
483 | // delete v0Lambda42; | |
484 | // delete v0Lambda24; | |
485 | // delete v0K0C; | |
486 | // delete v0GammaC; | |
487 | // delete v0Lambda42C; | |
488 | // delete v0Lambda24C; | |
489 | // } | |
490 | // ProcessPI0(); | |
491 | // } | |
492 | ||
493 | ||
494 | ||
495 | ||
496 | void AliTPCcalibV0::ProcessV0(Int_t ftype){ | |
497 | // | |
498 | // | |
499 | ||
500 | if (! fGammas) fGammas = new TObjArray(10); | |
501 | fGammas->Clear(); | |
502 | Int_t nV0s = fV0s->GetEntries(); | |
503 | if (nV0s==0) return; | |
504 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
505 | // | |
506 | for (Int_t ivertex=0; ivertex<nV0s; ivertex++){ | |
507 | AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex); | |
508 | AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track | |
509 | AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track | |
510 | ||
511 | const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam(); | |
512 | const AliExternalTrackParam * paramInPos = trackP->GetInnerParam(); | |
513 | ||
514 | if (!paramInPos) continue; // in case the inner paramters do not exist | |
515 | if (!paramInNeg) continue; | |
516 | // | |
517 | // | |
518 | // | |
519 | AliKFParticle *v0K0 = Fit(primVtx,v0,-211,211); | |
520 | AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11); | |
521 | AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211); | |
522 | AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212); | |
523 | //Set production vertex | |
524 | v0K0->SetProductionVertex( primVtx ); | |
525 | v0Gamma->SetProductionVertex( primVtx ); | |
526 | v0Lambda42->SetProductionVertex( primVtx ); | |
527 | v0Lambda24->SetProductionVertex( primVtx ); | |
528 | Double_t massK0, massGamma, massLambda42,massLambda24, massSigma; | |
529 | v0K0->GetMass( massK0,massSigma); | |
530 | v0Gamma->GetMass( massGamma,massSigma); | |
531 | v0Lambda42->GetMass( massLambda42,massSigma); | |
532 | v0Lambda24->GetMass( massLambda24,massSigma); | |
533 | Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF(); | |
534 | Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF(); | |
535 | Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF(); | |
536 | Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF(); | |
537 | // | |
538 | // Mass Contrained params | |
539 | // | |
540 | AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211); | |
541 | AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11); | |
542 | AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar | |
543 | AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda | |
544 | // | |
545 | v0K0C->SetProductionVertex( primVtx ); | |
546 | v0GammaC->SetProductionVertex( primVtx ); | |
547 | v0Lambda42C->SetProductionVertex( primVtx ); | |
548 | v0Lambda24C->SetProductionVertex( primVtx ); | |
549 | ||
550 | v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass()); | |
551 | v0GammaC->SetMassConstraint(0); | |
552 | v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass()); | |
553 | v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass()); | |
554 | // | |
555 | Double_t timeK0, sigmaTimeK0; | |
556 | Double_t timeLambda42, sigmaTimeLambda42; | |
557 | Double_t timeLambda24, sigmaTimeLambda24; | |
558 | v0K0C->GetLifeTime(timeK0, sigmaTimeK0); | |
559 | //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0); | |
560 | v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42); | |
561 | v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24); | |
562 | ||
563 | ||
564 | // | |
565 | Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF(); | |
566 | if (chi2K0C<0) chi2K0C=100; | |
567 | Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF(); | |
568 | if (chi2GammaC<0) chi2GammaC=100; | |
569 | Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF(); | |
570 | if (chi2Lambda42C<0) chi2Lambda42C=100; | |
571 | Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF(); | |
572 | if (chi2Lambda24C<0) chi2Lambda24C=100; | |
573 | // | |
574 | Float_t minChi2C=99; | |
575 | Int_t type =-1; | |
576 | if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;} | |
577 | if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;} | |
578 | if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;} | |
579 | if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;} | |
580 | Float_t minChi2=99; | |
581 | Int_t type0 =-1; | |
582 | if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;} | |
583 | if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;} | |
584 | if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;} | |
585 | if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;} | |
586 | ||
587 | // 0 is negative particle; 1 is positive particle | |
588 | Float_t betaGamma0 = 0; | |
589 | Float_t betaGamma1 = 0; | |
590 | ||
591 | switch (type) { | |
592 | case 0: | |
593 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); | |
594 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); | |
595 | break; | |
596 | case 1: | |
597 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass(); | |
598 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass(); | |
599 | break; | |
600 | case 2: | |
601 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass(); | |
602 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass(); | |
603 | break; | |
604 | case 3: | |
605 | betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass(); | |
606 | betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass(); | |
607 | break; | |
608 | } | |
609 | ||
610 | // cuts and histogram filling | |
611 | Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2 | |
612 | ||
613 | if (minChi2C < 2 && ftype == 1) { | |
614 | // | |
615 | if (chi2K0C < 10*minChi2C) numCand++; | |
616 | if (chi2GammaC < 10*minChi2C) numCand++; | |
617 | if (chi2Lambda42C < 10*minChi2C) numCand++; | |
618 | if (chi2Lambda24C < 10*minChi2C) numCand++; | |
619 | // | |
620 | if (numCand < 2) { | |
621 | if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal()); | |
622 | if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal()); | |
623 | } | |
624 | } | |
625 | ||
626 | // | |
627 | // | |
628 | // write output tree | |
629 | if (minChi2>50) continue; | |
630 | (*fDebugStream)<<"V0"<< | |
631 | "ftype="<<ftype<< | |
632 | "v0.="<<v0<< | |
633 | "trackN.="<<trackN<< | |
634 | "trackP.="<<trackP<< | |
635 | // | |
636 | "betaGamma0="<<betaGamma0<< | |
637 | "betaGamma1="<<betaGamma1<< | |
638 | // | |
639 | "type="<<type<< | |
640 | "chi2C="<<minChi2C<< | |
641 | "v0K0.="<<v0K0<< | |
642 | "v0Gamma.="<<v0Gamma<< | |
643 | "v0Lambda42.="<<v0Lambda42<< | |
644 | "v0Lambda24.="<<v0Lambda24<< | |
645 | // | |
646 | "chi20K0.="<<chi2K0<< | |
647 | "chi2Gamma.="<<chi2Gamma<< | |
648 | "chi2Lambda42.="<<chi2Lambda42<< | |
649 | "chi2Lambda24.="<<chi2Lambda24<< | |
650 | // | |
651 | "chi20K0c.="<<chi2K0C<< | |
652 | "chi2Gammac.="<<chi2GammaC<< | |
653 | "chi2Lambda42c.="<<chi2Lambda42C<< | |
654 | "chi2Lambda24c.="<<chi2Lambda24C<< | |
655 | // | |
656 | "v0K0C.="<<v0K0C<< | |
657 | "v0GammaC.="<<v0GammaC<< | |
658 | "v0Lambda42C.="<<v0Lambda42C<< | |
659 | "v0Lambda24C.="<<v0Lambda24C<< | |
660 | // | |
661 | "massK0="<<massK0<< | |
662 | "massGamma="<<massGamma<< | |
663 | "massLambda42="<<massLambda42<< | |
664 | "massLambda24="<<massLambda24<< | |
665 | // | |
666 | "timeK0="<<timeK0<< | |
667 | "timeLambda42="<<timeLambda42<< | |
668 | "timeLambda24="<<timeLambda24<< | |
669 | "\n"; | |
670 | if (type==1) fGammas->AddLast(v0); | |
671 | // | |
672 | // | |
673 | // | |
674 | delete v0K0; | |
675 | delete v0Gamma; | |
676 | delete v0Lambda42; | |
677 | delete v0Lambda24; | |
678 | delete v0K0C; | |
679 | delete v0GammaC; | |
680 | delete v0Lambda42C; | |
681 | delete v0Lambda24C; | |
682 | } | |
683 | ProcessPI0(); | |
684 | } | |
685 | ||
686 | ||
687 | ||
688 | void AliTPCcalibV0::ProcessPI0(){ | |
689 | // | |
690 | // | |
691 | // | |
692 | Int_t nentries = fGammas->GetEntries(); | |
693 | if (nentries<2) return; | |
694 | // | |
695 | Double_t m0[3], m1[3]; | |
696 | AliKFVertex primVtx(*(fESD->GetPrimaryVertex())); | |
697 | for (Int_t i0=0; i0<nentries; i0++){ | |
698 | AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0); | |
699 | v00->GetPxPyPz (m0[0], m0[1], m0[2]); | |
700 | AliKFParticle *p00 = Fit(primVtx, v00, 11,-11); | |
701 | p00->SetProductionVertex( primVtx ); | |
702 | p00->SetMassConstraint(0); | |
703 | // | |
704 | for (Int_t i1=i0; i1<nentries; i1++){ | |
705 | AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1); | |
706 | v01->GetPxPyPz (m1[0], m1[1], m1[2]); | |
707 | AliKFParticle *p01 = Fit(primVtx, v01, 11,-11); | |
708 | p01->SetProductionVertex( primVtx ); | |
709 | p01->SetMassConstraint(0); | |
710 | if (v00->GetIndex(0) != v01->GetIndex(0) && | |
711 | v00->GetIndex(1) != v01->GetIndex(1)){ | |
712 | AliKFParticle pi0( *p00,*p01); | |
713 | pi0.SetProductionVertex(primVtx); | |
714 | Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]); | |
715 | Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]); | |
716 | Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2]))); | |
717 | (*fDebugStream)<<"PI0"<< | |
718 | "v00.="<<v00<< | |
719 | "v01.="<<v01<< | |
720 | "mass="<<mass<< | |
721 | "p00.="<<p00<< | |
722 | "p01.="<<p01<< | |
723 | "pi0.="<<&pi0<< | |
724 | "\n"; | |
725 | } | |
726 | delete p01; | |
727 | } | |
728 | delete p00; | |
729 | } | |
730 | } | |
731 | ||
732 | ||
733 | ||
734 | ||
735 | ||
736 | AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & primVtx, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){ | |
737 | // | |
738 | // Make KF Particle | |
739 | // | |
740 | AliKFParticle p1( *(v0->GetParamN()), PDG1 ); | |
741 | AliKFParticle p2( *(v0->GetParamP()), PDG2 ); | |
742 | AliKFParticle *V0 = new AliKFParticle; | |
743 | Double_t x, y, z; | |
744 | v0->GetXYZ(x,y,z ); | |
745 | V0->SetVtxGuess(x,y,z); | |
746 | *(V0)+=p1; | |
747 | *(V0)+=p2; | |
748 | return V0; | |
749 | } | |
750 | ||
751 | ||
752 | ||
753 | Float_t AliTPCcalibV0::TPCBetheBloch(Float_t bg) | |
754 | { | |
755 | // | |
756 | // Bethe-Bloch energy loss formula | |
757 | // | |
758 | const Double_t kp1=0.76176e-1; | |
759 | const Double_t kp2=10.632; | |
760 | const Double_t kp3=0.13279e-4; | |
761 | const Double_t kp4=1.8631; | |
762 | const Double_t kp5=1.9479; | |
763 | Double_t dbg = (Double_t) bg; | |
764 | Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg); | |
765 | Double_t aa = TMath::Power(beta,kp4); | |
766 | Double_t bb = TMath::Power(1./dbg,kp5); | |
767 | bb=TMath::Log(kp3+bb); | |
768 | return ((Float_t)((kp2-aa-bb)*kp1/aa)); | |
769 | } | |
770 | ||
771 | ||
772 | ||
773 | ||
774 | ||
775 | ||
776 | ||
777 | TH2F * AliTPCcalibV0::GetHistograms() { | |
778 | // | |
779 | // | |
780 | // | |
781 | return fTPCdEdx; | |
782 | } | |
783 | ||
784 | ||
785 | ||
786 | ||
787 | void AliTPCcalibV0::BinLogX(TH2F *h) { | |
788 | // | |
789 | // | |
790 | // | |
791 | TAxis *axis = h->GetXaxis(); | |
792 | int bins = axis->GetNbins(); | |
793 | ||
794 | Double_t from = axis->GetXmin(); | |
795 | Double_t to = axis->GetXmax(); | |
796 | Double_t *new_bins = new Double_t[bins + 1]; | |
797 | ||
798 | new_bins[0] = from; | |
799 | Double_t factor = pow(to/from, 1./bins); | |
800 | ||
801 | for (int i = 1; i <= bins; i++) { | |
802 | new_bins[i] = factor * new_bins[i-1]; | |
803 | } | |
804 | axis->Set(bins, new_bins); | |
805 | delete new_bins; | |
806 | ||
807 | } | |
808 | ||
809 | ||
810 | ||
811 |