]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCcalibV0.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Corrected overlaps
[u/mrichter/AliRoot.git] / TPC / AliTPCcalibV0.cxx
CommitLineData
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
45ClassImp(AliTPCcalibV0)
46
47
48AliTPCcalibV0::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
79AliTPCcalibV0::~AliTPCcalibV0(){
80 //
81 //
82 //
83 delete fDebugStream;
84}
85
86
87
88void AliTPCcalibV0::ProofSlaveBegin(TList * output)
89{
90 // Called on PROOF - fill output list
91}
92
93
94void 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
109void 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
267void 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
497void 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
689void 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
737AliKFParticle * 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
754Float_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
778TH2F * AliTPCcalibV0::GetHistograms() {
779 //
780 //
781 //
782 return fTPCdEdx;
783}
784
785
786
787
788void 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