TPC/AliTPCcalibV0.cxx

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