]> 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
Coverty fixes in AliTPCCorrection
[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"
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
46ClassImp(AliTPCcalibV0)
47
48
49AliTPCcalibV0::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
74AliTPCcalibV0::~AliTPCcalibV0(){
75 //
76 //
77 //
10757ee9 78}
79
80
81
10757ee9 82
83
57e4988a 84void 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
99void 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
262void 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
492void 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
687void 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 738AliKFParticle * 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 753TH2F * AliTPCcalibV0::GetHistograms() {
754 //
755 //
756 //
757 return fTPCdEdx;
758}
759
760
761
762
763void 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);
4ce766eb 781 delete [] new_bins;
10757ee9 782
783}
784
785
786
787