PWG2/KINK/AliAnalysisKinkESDMC.cxx

   1 /**************************************************************************
   2  * Authors: Martha Spyropoulou-Stassinaki and the  members
   3  * of the Greek group at Physics Department of Athens University
   4  * Paraskevi Ganoti, Anastasia Belogianni and Filimon Roukoutakis
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------
  17 //                 AliAnalysisKinkESDMC class
  18 //       Example of an analysis task for kink topology study
  19 //      Kaons from kink topology are 'identified' in this code
  20 //-----------------------------------------------------------------
  21
  22 #include "TChain.h"
  23 #include "TTree.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TH1D.h"
  27 #include "TH2D.h"
  28 #include "TParticle.h"
  29 #include <TVector3.h>
  30 #include "TF1.h"
  31
  32 #include "AliAnalysisTaskSE.h"
  33 #include "AliAnalysisManager.h"
  34
  35 #include "AliESDEvent.h"
  36 #include "AliESDInputHandler.h"
  37
  38 #include "AliMCEventHandler.h"
  39 #include "AliMCEvent.h"
  40
  41 #include "AliAnalysisKinkESDMC.h"
  42
  43 #include "AliStack.h"
  44 #include "AliESDpid.h"
  45 #include "AliESDkink.h"
  46
  47 ClassImp(AliAnalysisKinkESDMC)
  48
  49 //________________________________________________________________________
  50 AliAnalysisKinkESDMC::AliAnalysisKinkESDMC()
  51   : AliAnalysisTaskSE(), fESD(0),fHistPtESD(),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
  52   , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
  53   fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0),  fgenPtEtR(0),fPtKink(0),  fptKink(0),
  54    fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0),
  55   fRpr(0),fZpr(0),
  56    fAngMomKC(0),
  57   fZvXv(0),fZvYv(0), fXvYv(0),  fPtPrKink(0), f1(0),  f2(0),
  58       fListOfHistos(0)
  59
  60 {
  61   // Constructor
  62
  63 }
  64
  65
  66 //________________________________________________________________________
  67 AliAnalysisKinkESDMC::AliAnalysisKinkESDMC(const char *name)
  68   : AliAnalysisTaskSE(name), fESD(0),fHistPtESD(0),fHistPt(0),fHistQtAll(0),fHistQt1(0),fHistQt2(0)
  69   , fHistPtKaon(0),fHistPtKPDG(0),fHistEta(0),fHistEtaK(0),fptKMC(0),fMultiplMC(0),fESDMult(0),fgenpt(0),frad(0),
  70   fKinkKaon(0), fKinkKaonBg(0), fM1kaon(0),  fgenPtEtR(0),fPtKink(0),  fptKink(0),
  71    fcodeH(0), fdcodeH(0), fAngMomK(0),fAngMomPi(0),
  72   fRpr(0),fZpr(0),
  73    fAngMomKC(0),
  74    fZvXv(0),fZvYv(0), fXvYv(0), fPtPrKink(0),  f1(0), f2(0),
  75       fListOfHistos(0)
  76   // , faod(0), fSignalTree(0)
  77
  78 {
  79   // Constructor
  80
  81   // Define input and output slots here
  82   // Input slot #0 works with a TChain
  83   DefineInput(0, TChain::Class());
  84   DefineOutput(1, TList::Class());
  85 }
  86
  87 //________________________________________________________________________
  88 void AliAnalysisKinkESDMC::ConnectInputData(Option_t *)
  89 {
  90   // Connect ESD or AOD here
  91   // Called once
  92
  93   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  94   if (!tree) {
  95     Printf("ERROR: Could not read chain from input slot 0");
  96   } else {
  97     // Disable all branches, we want to process only MC
  98     tree->SetBranchStatus("*", kTRUE);
  99    // tree->SetBranchStatus("fTracks.*", kTRUE);
 100    tree->SetBranchStatus("*Calo*", kFALSE);
 101    tree->SetBranchStatus("FMD", kFALSE);
 102
 103     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 104
 105     if (!esdH) {
 106       Printf("ERROR: Could not get ESDInputHandler");
 107     } else
 108       fESD = esdH->GetEvent();
 109   }
 110 }
 111
 112 //________________________________________________________________________
 113 void AliAnalysisKinkESDMC::CreateOutputObjects()
 114 {
 115   // Create histograms
 116   // Called once
 117
 118    f1=new TF1("f1","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",1.1,10.0);
 119    f1->SetParameter(0,0.493677);
 120    f1->SetParameter(1,0.9127037);
 121    f1->SetParameter(2,TMath::Pi());
 122
 123
 124    f2=new TF1("f2","((atan([0]*[1]*(1.0/(sqrt((x^2)*(1.0-([1]^2))-([0]^2)*([1]^2))))))*180.)/[2]",0.1,10.0);
 125    f2->SetParameter(0,0.13957018);
 126    f2->SetParameter(1,0.2731374);
 127    f2->SetParameter(2,TMath::Pi());
 128    //Open file  1= CAF
 129     //OpenFile(1);
 130
 131   fHistPtESD = new TH1F("fHistPtESD", "P_{T} distribution", 60, 0.0, 6.0);
 132   fHistPtESD->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 133   fHistPtESD->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 134   fHistPtESD->SetMarkerStyle(kFullCircle);
 135   fHistPt = new TH1F("fHistPt", "P_{T} distribution", 60, 0.0, 6.0);
 136   fHistQtAll = new TH1F("fHistQtAll", "Q_{T} distr All Kinks ",100, 0.0,.250);
 137   fHistQt1= new TH1F("fHistQt1", "Q_{T} distribution",100, 0.0,.300);
 138   fHistQt2= new TH1F("fHistQt2", "Q_{T} distribution",100, 0.0,.300);
 139   fHistPtKaon = new TH1F("fHistPtKaon", "P_{T}Kaon distribution", 60, 0.0, 6.0);
 140   fHistPtKPDG = new TH1F("fHistPtKPDG", "P_{T}Kaon distribution", 60, 0.0, 6.0);
 141   fHistEta= new TH1F("fHistEta", "Eta distribution", 26,-1.3, 1.3);
 142   fHistEtaK= new TH1F("fHistEtaK", "EtaK distribution", 26,-1.3, 1.3);
 143   fptKMC= new TH1F("fptKMC", "P_{T}Kaon generated", 60, 0.0, 6.0);
 144   fMultiplMC= new TH1F("fMultiplMC", "charge multiplicity MC", 60, 0.5,120.5);
 145   fESDMult= new TH1F("fESDMult", "charge multipliESD", 60, 0.5,120.5);
 146   fgenpt= new TH1F("fgenpt", "genpt   K distribution", 60, 0.0, 6.0);
 147   frad= new TH1F("frad", "radius  K generated",100,0.,400.);
 148   fKinkKaon= new TH1F("fKinkKaon", "P_{T}Kaon kinks identi", 60, 0.0, 6.0);
 149   fKinkKaonBg= new TH1F("fKinkKaonBg", "P_{T}Kaon kinks backgr", 60, 0.0, 6.0);
 150   fM1kaon= new TH1F("fM1kaon","Invar m(kaon) from kink->mu+netrino decay",80,0.0, 0.8);
 151   fgenPtEtR= new TH1F("fgenPtEtR", "P_{T}Kaon distribution", 60, 0.0, 6.0);
 152   fPtKink= new TH1F("fPtKink", "P_{T}Kaon Kink  bution", 60, 0.0, 6.0);
 153   fptKink= new TH1F("fptKink", "P_{T}Kaon Kink  bution", 60, 0.0, 6.0);
 154   fcodeH   = new TH2F("fcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 155   fdcodeH = new TH2F("fdcodeH", "code vrs dcode dist. kinks,K",100,0.,2500.,100,0.,2500.);
 156   fAngMomK= new TH2F("fAngMomK","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
 157   fAngMomPi= new TH2F("fAngMomPi","Decay angle vrs Mother Mom,Pi",100,0.0,5.0,80,0.,80.);
 158
 159   fRpr = new TH1D("fRpr", "rad distribution  PID pr",50,0.0, 2.5);
 160   fZpr = new TH1D("fZpr", "z distribution PID pr  ",60,-15.,15.);
 161   fAngMomKC= new TH2F("fAngMomKC","Decay angle vrs Mother Mom,K",100,0.0,5.0,80,0.,80.);
 162   fZvXv= new TH2F("fZvXv","Xv-Zv main vtx",60,-0.5,0.5,60, -15., 15.0);
 163   fZvYv= new TH2F("fZvYv","Yv-Zv main vtx",60,-0.5,0.5, 60, -15., 15.);
 164   fXvYv= new TH2F("fXvYv","Xv-Yv main vtx", 60,-1.5,1.5, 60, -1.5, 1.5);
 165   fPtPrKink=new TH1F("fPtPrKink","pt of ESD  kaonKink tracks", 60, 0.0, 6.0);
 166
 167    fListOfHistos=new TList();
 168
 169    fListOfHistos->Add(fHistPtESD);
 170    fListOfHistos->Add(fHistPt);
 171    fListOfHistos->Add(fHistQtAll);
 172    fListOfHistos->Add(fHistQt1);
 173    fListOfHistos->Add(fHistQt2);
 174    fListOfHistos->Add(fHistPtKaon);
 175    fListOfHistos->Add(fHistPtKPDG);
 176    fListOfHistos->Add(fHistEta);
 177    fListOfHistos->Add(fHistEtaK);
 178    fListOfHistos->Add(fptKMC);
 179    fListOfHistos->Add(fMultiplMC);
 180    fListOfHistos->Add(fESDMult);
 181    fListOfHistos->Add(fgenpt);
 182    fListOfHistos->Add(frad);
 183    fListOfHistos->Add(fKinkKaon);
 184    fListOfHistos->Add(fKinkKaonBg);
 185    fListOfHistos->Add(fM1kaon);
 186    fListOfHistos->Add(fgenPtEtR);
 187    fListOfHistos->Add(fPtKink);
 188    fListOfHistos->Add(fptKink);
 189    fListOfHistos->Add(fcodeH);
 190    fListOfHistos->Add(fdcodeH);
 191    fListOfHistos->Add(fdcodeH);
 192    fListOfHistos->Add(fAngMomK);
 193    fListOfHistos->Add(fAngMomPi);
 194    fListOfHistos->Add(fRpr);
 195    fListOfHistos->Add(fZpr);
 196    fListOfHistos->Add(fAngMomKC);
 197    fListOfHistos->Add(fZvXv);
 198    fListOfHistos->Add(fZvYv);
 199    fListOfHistos->Add(fXvYv);
 200    fListOfHistos->Add(fPtPrKink);
 201 }
 202
 203 //________________________________________________________________________
 204 void AliAnalysisKinkESDMC::Exec(Option_t *)
 205 {
 206   // Main loop
 207   // Called for each event
 208
 209   // Process MC truth, therefore we receive the AliAnalysisManager and ask it for the AliMCEventHandler
 210   // This handler can return the current MC event
 211
 212   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 213   if (!eventHandler) {
 214     Printf("ERROR: Could not retrieve MC event handler");
 215     return;
 216   }
 217
 218   AliMCEvent* mcEvent = eventHandler->MCEvent();
 219   if (!mcEvent) {
 220      Printf("ERROR: Could not retrieve MC event");
 221      return;
 222   }
 223
 224   Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 225
 226   if (!fESD) {
 227     Printf("ERROR: fESD not available");
 228     return;
 229   }
 230   AliStack* stack=mcEvent->Stack();
 231
 232 //primary tracks  in MC
 233          Int_t  nPrim = stack->GetNprimary();
 234 //
 235   const AliESDVertex *vertex=GetEventVertex(fESD);
 236     if(!vertex) return;
 237 //
 238   Double_t vpos[3];
 239   vertex->GetXYZ(vpos);
 240     fZpr->Fill(vpos[2]);
 241  fZvXv->Fill( vpos[1], vpos[2]);
 242  fZvYv->Fill( vpos[0], vpos[2]);
 243
 244   Double_t vtrack[3], ptrack[3];
 245
 246
 247
 248    Int_t nGoodTracks =  fESD->GetNumberOfTracks();
 249     fESDMult->Fill(nGoodTracks);
 250 //
 251 // track loop
 252   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 253     AliESDtrack* track = fESD->GetTrack(iTracks);
 254     if (!track) {
 255       Printf("ERROR: Could not receive track %d", iTracks);
 256       continue;
 257     }
 258
 259
 260     fHistPt->Fill(track->Pt());
 261
 262       //   pt cut at 300 MeV
 263         if ( (track->Pt())<.300)continue;
 264
 265     UInt_t status=track->GetStatus();
 266
 267   //  if((status&AliESDtrack::kITSrefit)==0) continue;
 268     if((status&AliESDtrack::kTPCrefit)==0) continue;
 269       if((track->GetTPCchi2()/track->GetTPCclusters(0))>3.5) continue;
 270
 271       Double_t extCovPos[15];
 272       track->GetExternalCovariance(extCovPos);
 273       if(extCovPos[0]>2) continue;
 274       if(extCovPos[2]>2) continue;
 275       if(extCovPos[5]>0.5) continue;
 276       if(extCovPos[9]>0.5) continue;
 277       if(extCovPos[14]>2) continue;
 278
 279
 280     track->GetXYZ(vtrack);
 281  fXvYv->Fill(vtrack[0],vtrack[1]);
 282
 283 // track momentum
 284      track->GetPxPyPz(ptrack);
 285
 286     TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
 287
 288     Double_t trackEta=trackMom.Eta();
 289
 290
 291
 292     Float_t nSigmaToVertex = GetSigmaToVertex(track);
 293
 294     Float_t bpos[2];
 295     Float_t bCovpos[3];
 296     track->GetImpactParameters(bpos,bCovpos);
 297
 298     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 299      Printf("Estimated b resolution lower or equal zero!");
 300      bCovpos[0]=0; bCovpos[2]=0;
 301     }
 302
 303     Float_t dcaToVertexXYpos = bpos[0];
 304     Float_t dcaToVertexZpos = bpos[1];
 305
 306     fRpr->Fill(dcaToVertexZpos);
 307
 308     if(nSigmaToVertex>=4) continue;
 309     if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
 310
 311
 312     Int_t label = track->GetLabel();
 313     label = TMath::Abs(label);
 314     TParticle * part = stack->Particle(label);
 315     if (!part) continue;
 316
 317  //  cut on eta
 318         if(  (TMath::Abs(trackEta )) > 1.1 ) continue;
 319     fHistPtESD->Fill(track->Pt());
 320
 321    // Add Kink analysis
 322
 323                 Int_t indexKinkPos=track->GetKinkIndex(0);
 324 //  loop on kinks
 325                 if(indexKinkPos<0){
 326                fPtKink->Fill(track->Pt()); /// pt from track
 327
 328         // select kink class
 329
 330           AliESDkink *kink=fESD->GetKink(TMath::Abs(indexKinkPos)-1);
 331 //
 332
 333           Int_t eSDfLabel1=kink->GetLabel(0);
 334           TParticle *particle1= stack->Particle(TMath::Abs(eSDfLabel1));
 335           Int_t code1= particle1->GetPdgCode();
 336
 337           Int_t eSDfLabeld=kink->GetLabel(1);
 338           TParticle *particled= stack->Particle(TMath::Abs(eSDfLabeld));
 339           Int_t dcode1= particled->GetPdgCode();
 340
 341            const TVector3 motherMfromKink(kink->GetMotherP());
 342            const TVector3 daughterMKink(kink->GetDaughterP());
 343            Float_t qT=kink->GetQt();
 344
 345            fHistEta->Fill(trackEta) ;  //   Eta distr
 346            fHistQtAll->Fill(qT) ;  //  Qt   distr
 347
 348            fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
 349
 350 //       fake kinks, small Qt and small kink angle
 351            if ( qT<0.012) fcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1) );
 352            if(qT<0.012) continue;  // remove fake kinks
 353
 354            Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 355
 356 //remove the double taracks
 357
 358            if( (kinkAngle<1.)  ) continue;
 359 //
 360            if( (qT>0.05)&& ( qT<0.25)  )   fHistQt2->Fill(qT);  // candidate kaon kinks
 361
 362            if (TMath::Abs(code1)==321)     fHistPtKPDG->Fill(track->Pt());  // ALL KAONS (pdg) inside kink sample
 363
 364 //          maximum decay angle at a given mother momentum
 365            Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
 366            Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
 367 //  two dimensional plot
 368                 if(TMath::Abs(code1)==321) fAngMomK->Fill(motherMfromKink.Mag(), kinkAngle);
 369                 if(TMath::Abs(code1)==211) fAngMomPi->Fill(motherMfromKink.Mag(), kinkAngle);
 370 //
 371 // invariant mass of mother track decaying to mu
 372          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
 373          Float_t p1XM= motherMfromKink.Px();
 374          Float_t p1YM= motherMfromKink.Py();
 375          Float_t p1ZM= motherMfromKink.Pz();
 376          Float_t p2XM= daughterMKink.Px();
 377          Float_t p2YM= daughterMKink.Py();
 378          Float_t p2ZM= daughterMKink.Pz();
 379          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 380          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 381
 382          fM1kaon->Fill(invariantMassKmu);
 383
 384 //  kaon selection from kinks
 385 if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackEta)<1.1)&&(invariantMassKmu<0.6)) {
 386                 fHistEtaK->Fill(trackEta);
 387
 388                 if(TMath::Abs(code1)==321) fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle);  // real kaons
 389
 390                 fHistQt1  ->Fill(qT) ;  //  Qt   distr
 391                 if ( (kinkAngle>maxDecAngKmu) && ( motherMfromKink.Mag()> 1.1 ) ) continue; // maximum angle selection
 392
 393                                     fHistPtKaon->Fill(track->Pt());   //all PID kink-kaon
 394
 395 // background inside the identified kaons, e.g KK  ,  Kp
 396          if((TMath::Abs(code1)==321)&&(( TMath::Abs(dcode1)) >=(TMath::Abs(code1)) )) fdcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1));
 397
 398          if((TMath::Abs(code1)==321) ) {
 399
 400          if ( label< nPrim ) fPtPrKink->Fill(track->Pt());
 401              fKinkKaon->Fill(track->Pt());        }
 402          else {
 403              fKinkKaonBg->Fill(track->Pt());
 404           }   // primary and secondary
 405
 406         }  //  kink selection
 407
 408
 409         }  //End Kink Information
 410
 411
 412   } //track loop
 413
 414   // loop over mc particles
 415
 416   // variable to count tracks
 417   Int_t nMCTracks = 0;
 418
 419 for (Int_t iMc = 0; iMc < nPrim; ++iMc)
 420   {
 421    // AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
 422
 423     TParticle* particle = stack->Particle(iMc);
 424
 425     if (!particle)
 426     {
 427     //  AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
 428       continue;
 429     }
 430
 431  //
 432     Float_t eta = particle->Eta();
 433
 434     if (TMath::Abs(eta) > 1.1)
 435       continue;
 436
 437             Double_t ptK = particle->Pt();
 438
 439    if( ptK <0.300) continue;
 440
 441       Float_t code = particle->GetPdgCode();
 442       if (code==321 || code==-321){
 443
 444               fptKMC->Fill(ptK);
 445
 446
 447             Int_t firstD=particle->GetFirstDaughter();
 448             Int_t lastD=particle->GetLastDaughter();
 449 //loop on secondaries
 450              for (Int_t k=firstD;k<lastD;k++) {
 451              TParticle*    daughter1=stack->Particle(k+1);
 452              Float_t dcode = daughter1->GetPdgCode();
 453
 454
 455
 456                  frad->Fill(daughter1->R());
 457
 458                  if (((daughter1->R())>110)&&((daughter1->R())<230) ){
 459         if (( ( code==321 )&& ( dcode ==-13  ))|| (( code == -321 )&& ( dcode == 13))) fgenPtEtR->Fill( ptK );//to muon
 460         if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211))) fgenPtEtR->Fill( ptK );//to pion
 461      //   if (( ( code==321 )&& ( dcode ==-11  ))|| (( code == -321 )&& ( dcode ==11))) fgenPtEtR->Fill( ptK );// to electr
 462
 463        }
 464 //  next only  to mu decay
 465                  if (((code==321)&&(dcode==-13))||((code==-321)&&(dcode==13))){
 466
 467
 468                     if (((daughter1->R())>110)&&((daughter1->R())<230)){
 469                     fgenpt->Fill(ptK);
 470                     }
 471                     }
 472        }
 473
 474       }
 475
 476
 477
 478     nMCTracks++;
 479   }// end of mc particle
 480
 481 //  printf("hello mc");
 482   fMultiplMC->Fill(nMCTracks);
 483
 484   PostData(1, fListOfHistos);
 485
 486 }
 487
 488 //________________________________________________________________________
 489 void AliAnalysisKinkESDMC::Terminate(Option_t *)
 490 {
 491   // Draw result to the screen
 492   // Called once at the end of the query
 493
 494 }
 495 //____________________________________________________________________//
 496
 497 Float_t AliAnalysisKinkESDMC::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 498   // Calculates the number of sigma to the vertex.
 499
 500   Float_t b[2];
 501   Float_t bRes[2];
 502   Float_t bCov[3];
 503
 504     esdTrack->GetImpactParameters(b,bCov);
 505
 506   if (bCov[0]<=0 || bCov[2]<=0) {
 507     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 508     bCov[0]=0; bCov[2]=0;
 509   }
 510   bRes[0] = TMath::Sqrt(bCov[0]);
 511   bRes[1] = TMath::Sqrt(bCov[2]);
 512
 513   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 514
 515   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 516
 517   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 518
 519   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 520
 521   return d;
 522   }
 523 //____________________________________________________________________//
 524
 525 const AliESDVertex* AliAnalysisKinkESDMC::GetEventVertex(const AliESDEvent* esd) const
 526   // Get the vertex from the ESD and returns it if the vertex is valid
 527
 528 {
 529   // Get the vertex
 530
 531   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 532
 533   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 534   else
 535   {
 536      vertex = esd->GetPrimaryVertexSPD();
 537      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 538      else
 539      return 0;
 540   }
 541 }