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 "AliESDEvent.h"
  27 #include "AliESDtrack.h"
  28 #include "AliESDfriend.h"
  29 #include "AliESDfriendTrack.h"
  30 #include "AliMathBase.h"
  31 #include "AliTPCseed.h"
  32 #include "AliTPCclusterMI.h"
  33
  34 #include "AliKFParticle.h"
  35 #include "AliKFVertex.h"
  36
  37 #include "AliTrackPointArray.h"
  38 #include "TCint.h"
  39 #include "AliTPCcalibV0.h"
  40 #include "AliV0.h"
  41
  42
  43
  44
  45
  46 ClassImp(AliTPCcalibV0)
  47
  48
  49 AliTPCcalibV0::AliTPCcalibV0() :
  50    AliTPCcalibBase(),
  51    fStack(0),
  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 {
  63   fPdg = new TDatabasePDG;
  64   // create output histograms
  65   fTPCdEdx   = new TH2F("TPCdEdX",  "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
  66   BinLogX(fTPCdEdx);
  67   fTPCdEdxPi = new TH2F("TPCdEdXPi","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
  68   fTPCdEdxEl = new TH2F("TPCdEdXEl","dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
  69   fTPCdEdxP  = new TH2F("TPCdEdXP", "dE/dX; BetaGamma; TPC signal (a.u.)", 1000, 0.1, 10000, 300, 0, 300);
  70
  71
  72 }
  73
  74 AliTPCcalibV0::~AliTPCcalibV0(){
  75   //
  76   //
  77   //
  78 }
  79
  80
  81
  82
  83
  84 void  AliTPCcalibV0::ProcessESD(AliESDEvent *esd, AliStack *stack){
  85   //
  86   //
  87   //
  88   fESD = esd;
  89   AliKFParticle::SetField(esd->GetMagneticField());
  90   MakeV0s();
  91   if (stack) {
  92     fStack = stack;
  93     MakeMC();
  94   }else{
  95     fStack =0;
  96   }
  97 }
  98
  99 void AliTPCcalibV0::MakeMC(){
 100   //
 101   // MC comparison
 102   //    1. Select interesting particles
 103   //    2. Assign the recosntructed particles
 104   //
 105   //1. Select interesting particles
 106   const Float_t kMinP   = 0.2;
 107   const Float_t kMinPt  = 0.1;
 108   const Float_t kMaxR   = 0.5;
 109   const Float_t kMaxTan = 1.2;
 110   const Float_t kMaxRad = 150;
 111   //
 112   if (!fParticles) fParticles = new TObjArray;
 113   TParticle *part=0;
 114   //
 115   Int_t entries = fStack->GetNtrack();
 116   for (Int_t ipart=0; ipart<entries; ipart++){
 117     part = fStack->Particle(ipart);
 118     if (!part) continue;
 119     if (part->P()<kMinP) continue;
 120     if (part->R()>kMaxR) continue;
 121     if (TMath::Abs(TMath::Tan(part->Theta()-TMath::Pi()*0.5))>kMaxTan) continue;
 122     Bool_t isInteresting = kFALSE;
 123     if (part->GetPdgCode()==22) isInteresting =kTRUE;
 124     if (part->GetPdgCode()==310) isInteresting =kTRUE;
 125     if (part->GetPdgCode()==111) isInteresting =kTRUE;
 126     if (TMath::Abs(part->GetPdgCode()==3122)) isInteresting =kTRUE;
 127
 128     //
 129     if (!isInteresting) continue;
 130     fParticles->AddLast(new TParticle(*part));
 131   }
 132   if (fParticles->GetEntries()<1) {
 133     return;
 134   }
 135   //
 136   //
 137   //
 138   Int_t sentries=fParticles->GetEntries();;
 139   for (Int_t ipart=0; ipart<sentries; ipart++){
 140     part = (TParticle*)fParticles->At(ipart);
 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       //
 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       }
 255     }
 256     fParticles->Delete();
 257   }
 258 }
 259
 260
 261
 262 void AliTPCcalibV0::MakeV0s(){
 263   //
 264   //
 265   //
 266   const Int_t kMinCluster=40;
 267   const Float_t kMinR    =0;
 268   if (! fV0s) fV0s = new TObjArray(10);
 269   fV0s->Clear();
 270   //
 271   // Old V0 finder
 272   //
 273   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
 274     AliESDv0 * v0 = fESD->GetV0(ivertex);
 275     if (v0->GetOnFlyStatus()) continue;
 276     fV0s->AddLast(v0);
 277   }
 278   ProcessV0(0);
 279   fV0s->Clear(0);
 280   //
 281   // MI V0 finder
 282   //
 283   for (Int_t ivertex=0; ivertex<fESD->GetNumberOfV0s(); ivertex++){
 284     AliESDv0 * v0 = fESD->GetV0(ivertex);
 285     if (!v0->GetOnFlyStatus()) continue;
 286     fV0s->AddLast(v0);
 287   }
 288   ProcessV0(1);
 289   fV0s->Clear();
 290   return;
 291   //
 292   // combinatorial
 293   //
 294   Int_t ntracks = fESD->GetNumberOfTracks();
 295   for (Int_t itrack0=0; itrack0<ntracks; itrack0++){
 296     AliESDtrack * track0 = fESD->GetTrack(itrack0);
 297     if (track0->GetSign()>0) continue;
 298     if ( track0->GetTPCNcls()<kMinCluster) continue;
 299     if (track0->GetKinkIndex(0)>0) continue;
 300     //
 301     for (Int_t itrack1=0; itrack1<ntracks; itrack1++){
 302       AliESDtrack * track1 = fESD->GetTrack(itrack1);
 303       if (track1->GetSign()<0) continue;
 304       if ( track1->GetTPCNcls()<kMinCluster) continue;
 305       if (track1->GetKinkIndex(0)>0) continue;
 306       //
 307       //      AliExternalTrackParam param0(*track0);
 308       // AliExternalTrackParam param1(*track1);
 309       AliV0 vertex;
 310       vertex.SetParamN(*track0);
 311       vertex.SetParamP(*track1);
 312       Float_t xyz[3];
 313       xyz[0] = fESD->GetPrimaryVertex()->GetXv();
 314       xyz[1] = fESD->GetPrimaryVertex()->GetYv();
 315       xyz[2] = fESD->GetPrimaryVertex()->GetZv();
 316       vertex.Update(xyz);
 317       if (vertex.GetRr()<kMinR) continue;
 318       if (vertex.GetDcaV0Daughters()>1.) continue;
 319       if (vertex.GetDcaV0Daughters()>0.3*vertex.GetRr()) continue;
 320       // if (vertex.GetPointAngle()<0.9) continue;
 321       vertex.SetIndex(0,itrack0);
 322       vertex.SetIndex(1,itrack1);
 323       fV0s->AddLast(new AliV0(vertex));
 324     }
 325   }
 326   ProcessV0(2);
 327   for (Int_t i=0;i<fV0s->GetEntries(); i++) delete fV0s->At(i);
 328   fV0s->Clear();
 329 }
 330
 331
 332
 333
 334
 335
 336 // void AliTPCcalibV0::ProcessV0(Int_t ftype){
 337 //   //
 338 //   //
 339 //   const Double_t ktimeK0     = 2.684;
 340 //   const Double_t ktimeLambda = 7.89;
 341
 342
 343 //   if (! fGammas) fGammas = new TObjArray(10);
 344 //   fGammas->Clear();
 345 //   Int_t nV0s  = fV0s->GetEntries();
 346 //   if (nV0s==0) return;
 347 //   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
 348 //   //
 349 //   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
 350 //     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
 351 //     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0));
 352 //     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1));
 353 //     //
 354 //     //
 355 //     //
 356 //     AliKFParticle *v0K0       = Fit(primVtx,v0,211,211);
 357 //     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
 358 //     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,2212,211);
 359 //     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,211,2212);
 360 //     //Set production vertex
 361 //     v0K0->SetProductionVertex( primVtx );
 362 //     v0Gamma->SetProductionVertex( primVtx );
 363 //     v0Lambda42->SetProductionVertex( primVtx );
 364 //     v0Lambda24->SetProductionVertex( primVtx );
 365 //     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
 366 //     v0K0->GetMass( massK0,massSigma);
 367 //     v0Gamma->GetMass( massGamma,massSigma);
 368 //     v0Lambda42->GetMass( massLambda42,massSigma);
 369 //     v0Lambda24->GetMass( massLambda24,massSigma);
 370 //     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
 371 //     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
 372 //     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
 373 //     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
 374 //     //
 375 //     // Mass Contrained params
 376 //     //
 377 //     AliKFParticle *v0K0C       = Fit(primVtx,v0,211,211);
 378 //     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
 379 //     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,2212,211);
 380 //     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,211,2212);
 381 //     //
 382 //     v0K0C->SetProductionVertex( primVtx );
 383 //     v0GammaC->SetProductionVertex( primVtx );
 384 //     v0Lambda42C->SetProductionVertex( primVtx );
 385 //     v0Lambda24C->SetProductionVertex( primVtx );
 386
 387 //     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
 388 //     v0GammaC->SetMassConstraint(0);
 389 //     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
 390 //     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
 391 //     //
 392 //     Double_t timeK0, sigmaTimeK0;
 393 //     Double_t timeLambda42, sigmaTimeLambda42;
 394 //     Double_t timeLambda24, sigmaTimeLambda24;
 395 //     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
 396 //     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
 397 //     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
 398 //     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
 399
 400
 401 //     //
 402 //     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
 403 //     if (chi2K0C<0) chi2K0C=100;
 404 //     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
 405 //     if (chi2GammaC<0) chi2GammaC=100;
 406 //     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
 407 //     if (chi2Lambda42C<0) chi2Lambda42C=100;
 408 //     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
 409 //     if (chi2Lambda24C<0) chi2Lambda24C=100;
 410 //     //
 411 //     Float_t  minChi2C=99;
 412 //     Int_t   type   =-1;
 413 //     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
 414 //     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
 415 //     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
 416 //     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
 417 //     Float_t  minChi2=99;
 418 //     Int_t   type0   =-1;
 419 //     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
 420 //     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
 421 //     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
 422 //     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
 423 //     Float_t betaGammaP  = trackN->GetP()/fPdg->GetParticle(-2212)->Mass();
 424 //     Float_t betaGammaPi = trackN->GetP()/fPdg->GetParticle(-211)->Mass();
 425 //     Float_t betaGammaEl = trackN->GetP()/fPdg->GetParticle(11)->Mass();
 426 //     Float_t dedxTeorP = BetheBlochAleph(betaGammaP);
 427 //     Float_t dedxTeorPi = BetheBlochAleph(betaGammaPi);;
 428 //     Float_t dedxTeorEl = BetheBlochAleph(betaGammaEl);;
 429 //     //
 430 //     //
 431 //     if (minChi2>50) continue;
 432 //     (*fDebugStream)<<"V0"<<
 433 //       "ftype="<<ftype<<
 434 //       "v0.="<<v0<<
 435 //       "trackN.="<<trackN<<
 436 //       "trackP.="<<trackP<<
 437 //       //
 438 //       "dedxTeorP="<<dedxTeorP<<
 439 //       "dedxTeorPi="<<dedxTeorPi<<
 440 //       "dedxTeorEl="<<dedxTeorEl<<
 441 //       //
 442 //       "type="<<type<<
 443 //       "chi2C="<<minChi2C<<
 444 //       "v0K0.="<<v0K0<<
 445 //       "v0Gamma.="<<v0Gamma<<
 446 //       "v0Lambda42.="<<v0Lambda42<<
 447 //       "v0Lambda24.="<<v0Lambda24<<
 448 //       //
 449 //       "chi20K0.="<<chi2K0<<
 450 //       "chi2Gamma.="<<chi2Gamma<<
 451 //       "chi2Lambda42.="<<chi2Lambda42<<
 452 //       "chi2Lambda24.="<<chi2Lambda24<<
 453 //       //
 454 //       "chi20K0c.="<<chi2K0C<<
 455 //       "chi2Gammac.="<<chi2GammaC<<
 456 //       "chi2Lambda42c.="<<chi2Lambda42C<<
 457 //       "chi2Lambda24c.="<<chi2Lambda24C<<
 458 //       //
 459 //       "v0K0C.="<<v0K0C<<
 460 //       "v0GammaC.="<<v0GammaC<<
 461 //       "v0Lambda42C.="<<v0Lambda42C<<
 462 //       "v0Lambda24C.="<<v0Lambda24C<<
 463 //       //
 464 //       "massK0="<<massK0<<
 465 //       "massGamma="<<massGamma<<
 466 //       "massLambda42="<<massLambda42<<
 467 //       "massLambda24="<<massLambda24<<
 468 //       //
 469 //       "timeK0="<<timeK0<<
 470 //       "timeLambda42="<<timeLambda42<<
 471 //       "timeLambda24="<<timeLambda24<<
 472 //       "\n";
 473 //     if (type==1) fGammas->AddLast(v0);
 474 //     //
 475 //     //
 476 //     //
 477 //     delete v0K0;
 478 //     delete v0Gamma;
 479 //     delete v0Lambda42;
 480 //     delete v0Lambda24;
 481 //     delete v0K0C;
 482 //     delete v0GammaC;
 483 //     delete v0Lambda42C;
 484 //     delete v0Lambda24C;
 485 //   }
 486 //   ProcessPI0();
 487 // }
 488
 489
 490
 491
 492 void AliTPCcalibV0::ProcessV0(Int_t ftype){
 493   //
 494   //
 495
 496   if (! fGammas) fGammas = new TObjArray(10);
 497   fGammas->Clear();
 498   Int_t nV0s  = fV0s->GetEntries();
 499   if (nV0s==0) return;
 500   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
 501   //
 502   for (Int_t ivertex=0; ivertex<nV0s; ivertex++){
 503     AliESDv0 * v0 = (AliESDv0*)fV0s->At(ivertex);
 504     AliESDtrack * trackN = fESD->GetTrack(v0->GetIndex(0)); // negative track
 505     AliESDtrack * trackP = fESD->GetTrack(v0->GetIndex(1)); // positive track
 506
 507     const AliExternalTrackParam * paramInNeg = trackN->GetInnerParam();
 508     const AliExternalTrackParam * paramInPos = trackP->GetInnerParam();
 509
 510     if (!paramInPos) continue; // in case the inner paramters do not exist
 511     if (!paramInNeg) continue;
 512     //
 513     //
 514     //
 515     AliKFParticle *v0K0       = Fit(primVtx,v0,-211,211);
 516     AliKFParticle *v0Gamma    = Fit(primVtx,v0,11,-11);
 517     AliKFParticle *v0Lambda42 = Fit(primVtx,v0,-2212,211);
 518     AliKFParticle *v0Lambda24 = Fit(primVtx,v0,-211,2212);
 519     //Set production vertex
 520     v0K0->SetProductionVertex( primVtx );
 521     v0Gamma->SetProductionVertex( primVtx );
 522     v0Lambda42->SetProductionVertex( primVtx );
 523     v0Lambda24->SetProductionVertex( primVtx );
 524     Double_t massK0, massGamma, massLambda42,massLambda24, massSigma;
 525     v0K0->GetMass( massK0,massSigma);
 526     v0Gamma->GetMass( massGamma,massSigma);
 527     v0Lambda42->GetMass( massLambda42,massSigma);
 528     v0Lambda24->GetMass( massLambda24,massSigma);
 529     Float_t chi2K0       = v0K0->GetChi2()/v0K0->GetNDF();
 530     Float_t chi2Gamma    = v0Gamma->GetChi2()/v0Gamma->GetNDF();
 531     Float_t chi2Lambda42 = v0Lambda42->GetChi2()/v0Lambda42->GetNDF();
 532     Float_t chi2Lambda24 = v0Lambda24->GetChi2()/v0Lambda24->GetNDF();
 533     //
 534     // Mass Contrained params
 535     //
 536     AliKFParticle *v0K0C       = Fit(primVtx,v0,-211,211);
 537     AliKFParticle *v0GammaC    = Fit(primVtx,v0,11,-11);
 538     AliKFParticle *v0Lambda42C = Fit(primVtx,v0,-2212,211); //lambdaBar
 539     AliKFParticle *v0Lambda24C = Fit(primVtx,v0,-211,2212); //lambda
 540     //
 541     v0K0C->SetProductionVertex( primVtx );
 542     v0GammaC->SetProductionVertex( primVtx );
 543     v0Lambda42C->SetProductionVertex( primVtx );
 544     v0Lambda24C->SetProductionVertex( primVtx );
 545
 546     v0K0C->SetMassConstraint(fPdg->GetParticle(310)->Mass());
 547     v0GammaC->SetMassConstraint(0);
 548     v0Lambda42C->SetMassConstraint(fPdg->GetParticle(-3122)->Mass());
 549     v0Lambda24C->SetMassConstraint(fPdg->GetParticle(3122)->Mass());
 550     //
 551     Double_t timeK0, sigmaTimeK0;
 552     Double_t timeLambda42, sigmaTimeLambda42;
 553     Double_t timeLambda24, sigmaTimeLambda24;
 554     v0K0C->GetLifeTime(timeK0, sigmaTimeK0);
 555     //v0K0Gamma->GetLifeTime(timeK0, sigmaTimeK0);
 556     v0Lambda42C->GetLifeTime(timeLambda42, sigmaTimeLambda42);
 557     v0Lambda24C->GetLifeTime(timeLambda24, sigmaTimeLambda24);
 558
 559
 560     //
 561     Float_t chi2K0C       = v0K0C->GetChi2()/v0K0C->GetNDF();
 562     if (chi2K0C<0) chi2K0C=100;
 563     Float_t chi2GammaC    = v0GammaC->GetChi2()/v0GammaC->GetNDF();
 564     if (chi2GammaC<0) chi2GammaC=100;
 565     Float_t chi2Lambda42C = v0Lambda42C->GetChi2()/v0Lambda42C->GetNDF();
 566     if (chi2Lambda42C<0) chi2Lambda42C=100;
 567     Float_t chi2Lambda24C = v0Lambda24C->GetChi2()/v0Lambda24C->GetNDF();
 568     if (chi2Lambda24C<0) chi2Lambda24C=100;
 569     //
 570     Float_t  minChi2C=99;
 571     Int_t   type   =-1;
 572     if (chi2K0C<minChi2C) { minChi2C= chi2K0C; type=0;}
 573     if (chi2GammaC<minChi2C) { minChi2C= chi2GammaC; type=1;}
 574     if (chi2Lambda42C<minChi2C) { minChi2C= chi2Lambda42C; type=2;}
 575     if (chi2Lambda24C<minChi2C) { minChi2C= chi2Lambda24C; type=3;}
 576     Float_t  minChi2=99;
 577     Int_t   type0   =-1;
 578     if (chi2K0<minChi2) { minChi2= chi2K0; type0=0;}
 579     if (chi2Gamma<minChi2) { minChi2= chi2Gamma; type0=1;}
 580     if (chi2Lambda42<minChi2) { minChi2= chi2Lambda42; type0=2;}
 581     if (chi2Lambda24<minChi2) { minChi2= chi2Lambda24; type0=3;}
 582
 583      // 0 is  negative particle; 1 is positive particle
 584     Float_t betaGamma0 = 0;
 585     Float_t betaGamma1 = 0;
 586
 587     switch (type) {
 588      case 0:
 589       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
 590       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
 591       break;
 592      case 1:
 593       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(11)->Mass();
 594       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(-11)->Mass();
 595       break;
 596      case 2:
 597       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-2212)->Mass();
 598       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(211)->Mass();
 599       break;
 600      case 3:
 601       betaGamma0 = paramInNeg->GetP()/fPdg->GetParticle(-211)->Mass();
 602       betaGamma1 = paramInPos->GetP()/fPdg->GetParticle(2212)->Mass();
 603       break;
 604     }
 605
 606     // cuts and histogram filling
 607     Int_t numCand = 0; // number of particle types which have a chi2 < 10*minChi2
 608
 609     if (minChi2C < 2 && ftype == 1) {
 610      //
 611      if (chi2K0C < 10*minChi2C) numCand++;
 612      if (chi2GammaC < 10*minChi2C) numCand++;
 613      if (chi2Lambda42C < 10*minChi2C) numCand++;
 614      if (chi2Lambda24C < 10*minChi2C) numCand++;
 615      //
 616      if (numCand < 2) {
 617       if (paramInNeg->GetP() > 0.4) fTPCdEdx->Fill(betaGamma0, trackN->GetTPCsignal());
 618       if (paramInPos->GetP() > 0.4) fTPCdEdx->Fill(betaGamma1, trackP->GetTPCsignal());
 619      }
 620     }
 621
 622     //
 623     //
 624     // write output tree
 625     if (minChi2>50) continue;
 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     }
 669     if (type==1) fGammas->AddLast(v0);
 670     //
 671     //
 672     //
 673     delete v0K0;
 674     delete v0Gamma;
 675     delete v0Lambda42;
 676     delete v0Lambda24;
 677     delete v0K0C;
 678     delete v0GammaC;
 679     delete v0Lambda42C;
 680     delete v0Lambda24C;
 681   }
 682   ProcessPI0();
 683 }
 684
 685
 686
 687 void AliTPCcalibV0::ProcessPI0(){
 688   //
 689   //
 690   //
 691   Int_t nentries = fGammas->GetEntries();
 692   if (nentries<2) return;
 693   //
 694   Double_t m0[3], m1[3];
 695   AliKFVertex primVtx(*(fESD->GetPrimaryVertex()));
 696   for (Int_t i0=0; i0<nentries; i0++){
 697     AliESDv0 *v00 = (AliESDv0*)fGammas->At(i0);
 698     v00->GetPxPyPz (m0[0], m0[1], m0[2]);
 699     AliKFParticle *p00 = Fit(primVtx, v00, 11,-11);
 700     p00->SetProductionVertex( primVtx );
 701     p00->SetMassConstraint(0);
 702     //
 703     for (Int_t i1=i0; i1<nentries; i1++){
 704       AliESDv0 *v01 = (AliESDv0*)fGammas->At(i1);
 705       v01->GetPxPyPz (m1[0], m1[1], m1[2]);
 706       AliKFParticle *p01 = Fit(primVtx, v01, 11,-11);
 707       p01->SetProductionVertex( primVtx );
 708       p01->SetMassConstraint(0);
 709       if (v00->GetIndex(0) != v01->GetIndex(0) &&
 710           v00->GetIndex(1) != v01->GetIndex(1)){
 711         AliKFParticle pi0( *p00,*p01);
 712         pi0.SetProductionVertex(primVtx);
 713         Double_t n1 = TMath::Sqrt (m0[0]*m0[0] + m0[1]*m0[1] + m0[2]*m0[2]);
 714         Double_t n2 = TMath::Sqrt (m1[0]*m1[0] + m1[1]*m1[1] + m1[2]*m1[2]);
 715         Double_t mass = TMath::Sqrt(2.*(n1*n2 - (m0[0]*m1[0] + m0[1]*m1[1] + m0[2]*m1[2])));
 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         }
 727       }
 728       delete p01;
 729     }
 730     delete p00;
 731   }
 732 }
 733
 734
 735
 736
 737
 738 AliKFParticle * AliTPCcalibV0::Fit(AliKFVertex & /*primVtx*/, AliESDv0 *v0, Int_t PDG1, Int_t PDG2){
 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
 753 TH2F * AliTPCcalibV0::GetHistograms() {
 754   //
 755   //
 756   //
 757  return fTPCdEdx;
 758 }
 759
 760
 761
 762
 763 void AliTPCcalibV0::BinLogX(TH2F *h) {
 764   //
 765   //
 766   //
 767    TAxis *axis = h->GetXaxis();
 768    int bins = axis->GetNbins();
 769
 770    Double_t from = axis->GetXmin();
 771    Double_t to = axis->GetXmax();
 772    Double_t *new_bins = new Double_t[bins + 1];
 773
 774    new_bins[0] = from;
 775    Double_t factor = pow(to/from, 1./bins);
 776
 777    for (int i = 1; i <= bins; i++) {
 778      new_bins[i] = factor * new_bins[i-1];
 779    }
 780    axis->Set(bins, new_bins);
 781    delete new_bins;
 782
 783 }
 784
 785
 786
 787