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