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