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