8 #include "TDatabasePDG.h"
9 #include "AliRunLoader.h"
16 #include "TLinearFitter.h"
18 #include "TTreeStream.h"
24 #include "AliTracker.h"
26 #include "AliESDtrack.h"
27 #include "AliESDfriend.h"
28 #include "AliESDfriendTrack.h"
29 #include "AliTPCseed.h"
30 #include "AliTPCclusterMI.h"
32 #include "AliKFParticle.h"
33 #include "AliKFVertex.h"
35 #include "AliTrackPointArray.h"
37 #include "AliTPCcalibV0.h"
44 ClassImp(AliTPCcalibV0)
47 AliTPCcalibV0::AliTPCcalibV0() :
63 G__SetCatchException(0);
64 fDebugStream = new TTreeSRedirector("V0debug.root");
65 fPdg = new TDatabasePDG;
68 // create output histograms
69 fTPCdEdx = new TH2F("TPCdEdX", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
71 fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
72 fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
73 fTPCdEdxP = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
78 AliTPCcalibV0::~AliTPCcalibV0(){
87 void AliTPCcalibV0::ProofSlaveBegin(TList * output)
89 // Called on PROOF - fill output list
93 void AliTPCcalibV0::ProcessESD(AliESD *esd, AliStack *stack){
98 AliKFParticle::SetField(esd->GetMagneticField());
108 void AliTPCcalibV0::MakeMC(){
111 // 1. Select interesting particles
112 // 2. Assign the recosntructed particles
114 //1. Select interesting particles
115 const Float_t kMinP = 0.2;
116 const Float_t kMinPt = 0.1;
117 const Float_t kMaxR = 0.5;
118 const Float_t kMaxTan = 1.2;
119 const Float_t kMaxRad = 150;
121 if (!fParticles) fParticles = new TObjArray;
124 Int_t entries = fStack->GetNtrack();
125 for (Int_t ipart=0; ipart<entries; ipart++){
126 part = fStack->Particle(ipart);
128 if (part->P()<kMinP) continue;
129 if (part->R()>kMaxR) continue;
130 if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue;
131 Bool_t isInteresting = kFALSE;
132 if (part->GetPdgCode()==22) isInteresting =kTRUE;
133 if (part->GetPdgCode()==310) isInteresting =kTRUE;
134 if (part->GetPdgCode()==111) isInteresting =kTRUE;
135 if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE;
138 if (!isInteresting) continue;
139 fParticles->AddLast(new TParticle(*part));
141 if (fParticles->GetEntries()<1) {
147 Int_t sentries=fParticles->GetEntries();;
148 for (Int_t ipart=0; ipart<sentries; ipart++){
149 TParticle *part = (TParticle*)fParticles->At(ipart);
155 Int_t id0 = part->GetDaughter(0);
156 Int_t id1 = part->GetDaughter(1);
157 if (id0>=fStack->GetNtrack() ) id0*=-1;
158 if (id1>=fStack->GetNtrack() ) id1*=-1;
159 Bool_t findable = kTRUE;
160 if (id0<0 || id1<0) findable = kFALSE;
163 p0 = fStack->Particle(id0);
164 if (p0->R()>kMaxRad) findable = kFALSE;
165 if (p0->Pt()<kMinPt) findable = kFALSE;
166 if (p0->Vz()>250) findable= kFALSE;
167 if (TMath::Abs(TMath::Tan(p0->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
168 if (fPdg->GetParticle(p0->GetPdgCode())==0) findable =kFALSE;
170 if (fPdg->GetParticle(p0->GetPdgCode())->Charge()==0) charge++;
172 p1 = fStack->Particle(id1);
173 if (p1->R()>kMaxRad) findable = kFALSE;
174 if (p1->Pt()<kMinPt) findable = kFALSE;
175 if (TMath::Abs(p1->Vz())>250) findable= kFALSE;
176 if (TMath::Abs(TMath::Tan(p1->Theta()-TMath::Pi()*0.5))>2) findable=kFALSE;
177 if (fPdg->GetParticle(p1->GetPdgCode())==0) findable = kFALSE;
179 if (fPdg->GetParticle(p1->GetPdgCode())->Charge()==0) charge++;
182 // (*fDebugStream)<<"MC0"<<
184 // "findable="<<findable<<
188 if (!findable) continue;
189 Float_t minpt = TMath::Min(p0->Pt(), p1->Pt());
194 AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
195 for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
196 AliESDv0 * v0 = fESD->GetV0(ivertex);
197 // select coresponding track
198 AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
199 if (TMath::Abs(trackN->GetLabel())!=id0 && TMath::Abs(trackN->GetLabel())!=id1) continue;
200 AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
201 if (TMath::Abs(trackP->GetLabel())!=id0 && TMath::Abs(trackP->GetLabel())!=id1) continue;
202 TParticle *pn = fStack->Particle(TMath::Abs(trackN->GetLabel()));
203 TParticle *pp = fStack->Particle(TMath::Abs(trackP->GetLabel()));
206 if ( v0->GetOnFlyStatus()) nnew++;
207 if (!v0->GetOnFlyStatus()) nold++;
208 if (part->GetPdgCode()==22 && TMath::Abs(pn->GetPdgCode())==11 && TMath::Abs(pp->GetPdgCode())==11)
210 if (part->GetPdgCode()==310 && TMath::Abs(pn->GetPdgCode())==211 && TMath::Abs(pp->GetPdgCode())==211)
212 if (part->GetPdgCode()==3122){
213 if (TMath::Abs(pn->GetPdgCode())==210 ) type=2;
216 AliKFParticle *v0kf = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
217 v0kf->SetProductionVertex( primVtx );
218 AliKFParticle *v0kfc = Fit(primVtx,v0,pn->GetPdgCode(),pp->GetPdgCode());
219 v0kfc->SetProductionVertex( primVtx );
220 v0kfc->SetMassConstraint(fPdg->GetParticle(part->GetPdgCode())->Mass());
221 Float_t chi2 = v0kf->GetChi2();
222 Float_t chi2C = v0kf->GetChi2();
225 (*fDebugStream)<<"MCRC"<<
247 (*fDebugStream)<<"MC"<<
260 fParticles->Delete();
266 void AliTPCcalibV0::MakeV0s(){
270 const Int_t kMinCluster=40;
271 const Float_t kMinR =0;
272 if (! fV0s) fV0s = new TObjArray(10);
277 for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
278 AliESDv0 * v0 = fESD->GetV0(ivertex);
279 if (v0->GetOnFlyStatus()) continue;
287 for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
288 AliESDv0 * v0 = fESD->GetV0(ivertex);
289 if (!v0->GetOnFlyStatus()) continue;
298 Int_t ntracks = fESD->GetNumberOfTracks();
299 for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
300 AliESDtrack * track0 = fESD->GetTrack(itrack0);
301 if (track0->GetSign()>0) continue;
302 if ( track0->GetTPCNcls()<kMinCluster) continue;
303 if (track0->GetKinkIndex(0)>0) continue;
305 for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
306 AliESDtrack * track1 = fESD->GetTrack(itrack1);
307 if (track1->GetSign()<0) continue;
308 if ( track1->GetTPCNcls()<kMinCluster) continue;
309 if (track1->GetKinkIndex(0)>0) continue;
311 // AliExternalTrackParam param0(*track0);
312 // AliExternalTrackParam param1(*track1);
314 vertex.SetParamN(*track0);
315 vertex.SetParamP(*track1);
317 xyz[0] = fESD->GetPrimaryVertex()->GetXv();
318 xyz[1] = fESD->GetPrimaryVertex()->GetYv();
319 xyz[2] = fESD->GetPrimaryVertex()->GetZv();
321 if (vertex.GetRr()<kMinR) continue;
322 if (vertex.GetDcaV0Daughters()>1.) continue;
323 if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue;
324 // if (vertex.GetPointAngle()<0.9) continue;
325 vertex.SetIndex(0,itrack0);
326 vertex.SetIndex(1,itrack1);
327 fV0s->AddLast(new AliV0(vertex));
331 for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i);
340 // void AliTPCcalibV0::ProcessV0(Int_t ftype){
343 // const Double_t ktimeK0 = 2.684;
344 // const Double_t ktimeLambda = 7.89;
347 // if (! fGammas) fGammas = new TObjArray(10);
349 // Int_t nV0s = fV0s->GetEntries();
350 // if (nV0s==0) return;
351 // AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
353 // for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
354 // AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
355 // AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
356 // AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
360 // AliKFParticle *v0K0 = Fit(primVtx,v0,211,211);
361 // AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11);
362 // AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211);
363 // AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212);
364 // //Set production vertex
365 // v0K0->SetProductionVertex( primVtx );
366 // v0Gamma->SetProductionVertex( primVtx );
367 // v0Lambda42->SetProductionVertex( primVtx );
368 // v0Lambda24->SetProductionVertex( primVtx );
369 // Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
370 // v0K0->GetMass( massK0,massSigma);
371 // v0Gamma->GetMass( massGamma,massSigma);
372 // v0Lambda42->GetMass( massLambda42,massSigma);
373 // v0Lambda24->GetMass( massLambda24,massSigma);
374 // Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF();
375 // Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF();
376 // Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
377 // Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
379 // // Mass Contrained params
381 // AliKFParticle *v0K0C = Fit(primVtx,v0,211,211);
382 // AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11);
383 // AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211);
384 // AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212);
386 // v0K0C->SetProductionVertex( primVtx );
387 // v0GammaC->SetProductionVertex( primVtx );
388 // v0Lambda42C->SetProductionVertex( primVtx );
389 // v0Lambda24C->SetProductionVertex( primVtx );
391 // v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
392 // v0GammaC->SetMassConstraint(0);
393 // v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
394 // v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
396 // Double_t timeK0, sigmaTimeK0;
397 // Double_t timeLambda42, sigmaTimeLambda42;
398 // Double_t timeLambda24, sigmaTimeLambda24;
399 // v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
400 // //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
401 // v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
402 // v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
406 // Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF();
407 // if (chi2K0C<0) chi2K0C=100;
408 // Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF();
409 // if (chi2GammaC<0) chi2GammaC=100;
410 // Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
411 // if (chi2Lambda42C<0) chi2Lambda42C=100;
412 // Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
413 // if (chi2Lambda24C<0) chi2Lambda24C=100;
415 // Float_t minChi2C=99;
417 // if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
418 // if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
419 // if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
420 // if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
421 // Float_t minChi2=99;
423 // if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
424 // if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
425 // if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
426 // if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
427 // Float_t betaGammaP = trackN->GetP()/fPdg->GetParticle(-2212)->Mass();
428 // Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass();
429 // Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass();
430 // Float_t dedxTeorP = TPCBetheBloch(betaGammaP);
431 // Float_t dedxTeorPi = TPCBetheBloch(betaGammaPi);;
432 // Float_t dedxTeorEl = TPCBetheBloch(betaGammaEl);;
435 // if (minChi2>50) continue;
436 // (*fDebugStream)<<"V0"<<
439 // "trackN.="<<trackN<<
440 // "trackP.="<<trackP<<
442 // "dedxTeorP="<<dedxTeorP<<
443 // "dedxTeorPi="<<dedxTeorPi<<
444 // "dedxTeorEl="<<dedxTeorEl<<
447 // "chi2C="<<minChi2C<<
449 // "v0Gamma.="<<v0Gamma<<
450 // "v0Lambda42.="<<v0Lambda42<<
451 // "v0Lambda24.="<<v0Lambda24<<
453 // "chi20K0.="<<chi2K0<<
454 // "chi2Gamma.="<<chi2Gamma<<
455 // "chi2Lambda42.="<<chi2Lambda42<<
456 // "chi2Lambda24.="<<chi2Lambda24<<
458 // "chi20K0c.="<<chi2K0C<<
459 // "chi2Gammac.="<<chi2GammaC<<
460 // "chi2Lambda42c.="<<chi2Lambda42C<<
461 // "chi2Lambda24c.="<<chi2Lambda24C<<
463 // "v0K0C.="<<v0K0C<<
464 // "v0GammaC.="<<v0GammaC<<
465 // "v0Lambda42C.="<<v0Lambda42C<<
466 // "v0Lambda24C.="<<v0Lambda24C<<
468 // "massK0="<<massK0<<
469 // "massGamma="<<massGamma<<
470 // "massLambda42="<<massLambda42<<
471 // "massLambda24="<<massLambda24<<
473 // "timeK0="<<timeK0<<
474 // "timeLambda42="<<timeLambda42<<
475 // "timeLambda24="<<timeLambda24<<
477 // if (type==1) fGammas->AddLast(v0);
483 // delete v0Lambda42;
484 // delete v0Lambda24;
487 // delete v0Lambda42C;
488 // delete v0Lambda24C;
496 void AliTPCcalibV0::ProcessV0(Int_t ftype){
500 if (! fGammas) fGammas = new TObjArray(10);
502 Int_t nV0s = fV0s->GetEntries();
504 AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
506 for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
507 AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
508 AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track
509 AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track
511 const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam();
512 const AliExternalTrackParam * paramInPos = trackP->GetInnerParam();
514 if (!paramInPos) continue; // in case the inner paramters do not exist
515 if (!paramInNeg) continue;
519 AliKFParticle *v0K0 = Fit(primVtx,v0,-211,211);
520 AliKFParticle *v0Gamma = Fit(primVtx,v0,11,-11);
521 AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211);
522 AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212);
523 //Set production vertex
524 v0K0->SetProductionVertex( primVtx );
525 v0Gamma->SetProductionVertex( primVtx );
526 v0Lambda42->SetProductionVertex( primVtx );
527 v0Lambda24->SetProductionVertex( primVtx );
528 Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
529 v0K0->GetMass( massK0,massSigma);
530 v0Gamma->GetMass( massGamma,massSigma);
531 v0Lambda42->GetMass( massLambda42,massSigma);
532 v0Lambda24->GetMass( massLambda24,massSigma);
533 Float_t chi2K0 = v0K0->GetChi2()/v0K0->GetNDF();
534 Float_t chi2Gamma = v0Gamma->GetChi2()/v0Gamma->GetNDF();
535 Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
536 Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
538 // Mass Contrained params
540 AliKFParticle *v0K0C = Fit(primVtx,v0,-211,211);
541 AliKFParticle *v0GammaC = Fit(primVtx,v0,11,-11);
542 AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar
543 AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda
545 v0K0C->SetProductionVertex( primVtx );
546 v0GammaC->SetProductionVertex( primVtx );
547 v0Lambda42C->SetProductionVertex( primVtx );
548 v0Lambda24C->SetProductionVertex( primVtx );
550 v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
551 v0GammaC->SetMassConstraint(0);
552 v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
553 v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
555 Double_t timeK0, sigmaTimeK0;
556 Double_t timeLambda42, sigmaTimeLambda42;
557 Double_t timeLambda24, sigmaTimeLambda24;
558 v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
559 //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
560 v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
561 v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
565 Float_t chi2K0C = v0K0C->GetChi2()/v0K0C->GetNDF();
566 if (chi2K0C<0) chi2K0C=100;
567 Float_t chi2GammaC = v0GammaC->GetChi2()/v0GammaC->GetNDF();
568 if (chi2GammaC<0) chi2GammaC=100;
569 Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
570 if (chi2Lambda42C<0) chi2Lambda42C=100;
571 Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
572 if (chi2Lambda24C<0) chi2Lambda24C=100;
576 if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
577 if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
578 if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
579 if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
582 if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
583 if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
584 if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
585 if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
587 // 0 is negative particle; 1 is positive particle
588 Float_t betaGamma0 = 0;
589 Float_t betaGamma1 = 0;
593 betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
594 betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
597 betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass();
598 betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass();
601 betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass();
602 betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
605 betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
606 betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass();
610 // cuts and histogram filling
611 Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2
613 if (minChi2C < 2 && ftype == 1) {
615 if (chi2K0C < 10*minChi2C) numCand++;
616 if (chi2GammaC < 10*minChi2C) numCand++;
617 if (chi2Lambda42C < 10*minChi2C) numCand++;
618 if (chi2Lambda24C < 10*minChi2C) numCand++;
621 if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal());
622 if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal());
629 if (minChi2>50) continue;
630 (*fDebugStream)<<"V0"<<
636 "betaGamma0="<<betaGamma0<<
637 "betaGamma1="<<betaGamma1<<
642 "v0Gamma.="<<v0Gamma<<
643 "v0Lambda42.="<<v0Lambda42<<
644 "v0Lambda24.="<<v0Lambda24<<
646 "chi20K0.="<<chi2K0<<
647 "chi2Gamma.="<<chi2Gamma<<
648 "chi2Lambda42.="<<chi2Lambda42<<
649 "chi2Lambda24.="<<chi2Lambda24<<
651 "chi20K0c.="<<chi2K0C<<
652 "chi2Gammac.="<<chi2GammaC<<
653 "chi2Lambda42c.="<<chi2Lambda42C<<
654 "chi2Lambda24c.="<<chi2Lambda24C<<
657 "v0GammaC.="<<v0GammaC<<
658 "v0Lambda42C.="<<v0Lambda42C<<
659 "v0Lambda24C.="<<v0Lambda24C<<
662 "massGamma="<<massGamma<<
663 "massLambda42="<<massLambda42<<
664 "massLambda24="<<massLambda24<<
667 "timeLambda42="<<timeLambda42<<
668 "timeLambda24="<<timeLambda24<<
670 if (type==1) fGammas->AddLast(v0);
688 void AliTPCcalibV0::ProcessPI0(){
692 Int_t nentries = fGammas->GetEntries();
693 if (nentries<2) return;
695 Double_t m0[3], m1[3];
696 AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
697 for (Int_t i0=0; i0<nentries; i0++){
698 AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0);
699 v00->GetPxPyPz (m0[0], m0[1], m0[2]);
700 AliKFParticle *p00 = Fit(primVtx, v00, 11,-11);
701 p00->SetProductionVertex( primVtx );
702 p00->SetMassConstraint(0);
704 for (Int_t i1=i0; i1<nentries; i1++){
705 AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1);
706 v01->GetPxPyPz (m1[0], m1[1], m1[2]);
707 AliKFParticle *p01 = Fit(primVtx, v01, 11,-11);
708 p01->SetProductionVertex( primVtx );
709 p01->SetMassConstraint(0);
710 if (v00->GetIndex(0) != v01->GetIndex(0) &&
711 v00->GetIndex(1) != v01->GetIndex(1)){
712 AliKFParticle pi0( *p00,*p01);
713 pi0.SetProductionVertex(primVtx);
714 Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]);
715 Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]);
716 Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2])));
717 (*fDebugStream)<<"PI0"<<
736 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & primVtx, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
740 AliKFParticle p1( *(v0->GetParamN()), PDG1 );
741 AliKFParticle p2( *(v0->GetParamP()), PDG2 );
742 AliKFParticle *V0 = new AliKFParticle;
745 V0->SetVtxGuess(x,y,z);
753 Float_t AliTPCcalibV0::TPCBetheBloch(Float_t bg)
756 // Bethe-Bloch energy loss formula
758 const Double_t kp1=0.76176e-1;
759 const Double_t kp2=10.632;
760 const Double_t kp3=0.13279e-4;
761 const Double_t kp4=1.8631;
762 const Double_t kp5=1.9479;
763 Double_t dbg = (Double_t) bg;
764 Double_t beta = dbg/TMath::Sqrt(1.+dbg*dbg);
765 Double_t aa = TMath::Power(beta,kp4);
766 Double_t bb = TMath::Power(1./dbg,kp5);
767 bb=TMath::Log(kp3+bb);
768 return ((Float_t)((kp2-aa-bb)*kp1/aa));
777 TH2F * AliTPCcalibV0::GetHistograms() {
787 void AliTPCcalibV0::BinLogX(TH2F *h) {
791 TAxis *axis = h->GetXaxis();
792 int bins = axis->GetNbins();
794 Double_t from = axis->GetXmin();
795 Double_t to = axis->GetXmax();
796 Double_t *new_bins = new Double_t[bins + 1];
799 Double_t factor = pow(to/from, 1./bins);
801 for (int i = 1; i <= bins; i++) {
802 new_bins[i] = factor * new_bins[i-1];
804 axis->Set(bins, new_bins);