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