Rulechecker-complying update from P.Ganoti (pganoti@phys.uoa.gr)
[u/mrichter/AliRoot.git] / 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
199   Double_t vtrack[3], ptrack[3];
200   
201      
202  // Int_t nESDTracK = 0;
203
204    Int_t nGoodTracks =  esd->GetNumberOfTracks();
205     fESDMult->Fill(nGoodTracks);
206      
207 //
208 // track loop
209    for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
210 // loop only on primary tracks
211
212     AliESDtrack* track = esd->GetTrack(iTracks);
213     if (!track) {
214       Printf("ERROR: Could not receive track %d", iTracks);
215       continue;
216     }
217     
218
219     fHistPt->Fill(track->Pt());
220
221     
222     Int_t label = track->GetLabel();
223     label = TMath::Abs(label);
224     TParticle * part = stack->Particle(label);
225     if (!part) continue;
226 // loop only on Primary tracks
227       if (label > nPrim) continue; /// primary tracks only   , 26/8/09  EFF study
228
229       //    pt cut at 300 MeV
230         if ( (track->Pt())<.300)continue;
231
232 //    UInt_t status=track->GetStatus();
233
234   //  if((status&AliESDtrack::kITSrefit)==0) continue;
235   //  if((status&AliESDtrack::kTPCrefit)==0) continue;
236   // ayksanei to ratio BG/real    if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue;
237      if((track->GetTPCchi2()/track->GetTPCclusters(0))>4.5) continue;
238
239       Double_t extCovPos[15];
240       track->GetExternalCovariance(extCovPos);    
241     //  if(extCovPos[0]>2) continue;
242     //  if(extCovPos[2]>2) continue;    
243     //  if(extCovPos[5]>0.5) continue;  
244     //  if(extCovPos[9]>0.5) continue;
245     //  if(extCovPos[14]>2) continue;
246
247
248     track->GetXYZ(vtrack);
249  fXvYv->Fill(vtrack[0],vtrack[1]);  
250  fZvYv->Fill(vtrack[0],vtrack[2]);  
251  fZvXv->Fill(vtrack[1],vtrack[2]);  
252
253 // track momentum
254      track->GetPxPyPz(ptrack);
255     
256     TVector3 trackMom(ptrack[0],ptrack[1],ptrack[2]);
257     
258     Double_t trackEta=trackMom.Eta();
259     
260     
261     
262     Float_t nSigmaToVertex = GetSigmaToVertex(track);  
263       
264     Float_t bpos[2];
265     Float_t bCovpos[3];
266     track->GetImpactParameters(bpos,bCovpos);
267     
268     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
269      Printf("Estimated b resolution lower or equal zero!");
270      bCovpos[0]=0; bCovpos[2]=0;
271     }
272
273     Float_t dcaToVertexXYpos = bpos[0];
274     Float_t dcaToVertexZpos = bpos[1];
275     
276     fRpr->Fill(dcaToVertexZpos);
277 //  test for secondary kinks   , 5/7/2009  
278     if(nSigmaToVertex>=4) continue;
279  //    if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;  //arxiko 
280      if((dcaToVertexXYpos>4.0)||(dcaToVertexZpos>5.0)) continue;   // 27/8
281     
282
283  //  cut on eta 
284         if(  (TMath::Abs(trackEta )) > 0.9 ) continue;
285     fHistPtESD->Fill(track->Pt());
286
287    // Add Kink analysis
288    
289                 Int_t indexKinkPos=track->GetKinkIndex(0);
290 //  loop on kinks
291                 if(indexKinkPos<0){
292                fPtKink->Fill(track->Pt()); /// pt from track
293
294         // select kink class    
295
296           AliESDkink *kink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
297 //
298         
299           Int_t eSDfLabel1=kink->GetLabel(0);
300           TParticle *particle1= stack->Particle(TMath::Abs(eSDfLabel1));
301           Int_t code1= particle1->GetPdgCode();
302           
303           Int_t eSDfLabeld=kink->GetLabel(1);
304           TParticle *particled= stack->Particle(TMath::Abs(eSDfLabeld));
305           Int_t dcode1= particled->GetPdgCode();
306
307            const TVector3 motherMfromKink(kink->GetMotherP());
308            const TVector3 daughterMKink(kink->GetDaughterP());
309            Float_t qT=kink->GetQt();
310
311            fHistQtAll->Fill(qT) ;  //  Qt   distr
312                   
313            fptKink->Fill(motherMfromKink.Pt()); /// pt from kink
314
315            Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
316
317 //       fake kinks, small Qt and small kink angle
318     if(( kinkAngle<1.)&&(TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13)) fHistQt1  ->Fill(qT) ;  //  Qt   distr
319
320            if(qT<0.012) continue;  // remove fake kinks
321
322 //remove the double taracks 
323            if( (kinkAngle<1.)  ) continue;
324
325 //
326
327       if((kink->GetR()>120.)&&(kink->GetR()<220.)&&(TMath::Abs(trackEta)<0.9)&&(label<nPrim)) {
328          if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||    
329            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))  ||    
330            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))  ) { 
331                fHistPtKPDG->Fill(track->Pt());  // ALL KAONS (pdg) inside ESD  kink sample
332            fHistEta->Fill(trackEta) ;  //   Eta distr of PDG kink ESD  kaons
333     fMultESDK->Fill(nGoodTracks);
334       fHistQt2->Fill(qT);  // PDG ESD kaons            
335      }
336      }
337
338 //          maximum decay angle at a given mother momentum
339            Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
340            Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
341 //  two dimensional plot 
342                 if(TMath::Abs(code1)==321) fAngMomK->Fill(motherMfromKink.Mag(), kinkAngle); 
343                 if(TMath::Abs(code1)==211) fAngMomPi->Fill(motherMfromKink.Mag(), kinkAngle); 
344 //
345 // invariant mass of mother track decaying to mu
346          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
347          Float_t p1XM= motherMfromKink.Px();
348          Float_t p1YM= motherMfromKink.Py();
349          Float_t p1ZM= motherMfromKink.Pz();
350          Float_t p2XM= daughterMKink.Px();
351          Float_t p2YM= daughterMKink.Py();
352          Float_t p2ZM= daughterMKink.Pz();
353          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
354          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
355   
356          fM1kaon->Fill(invariantMassKmu);
357
358 //  kaon selection from kinks
359 //if((kinkAngle>maxDecAngpimu)&&(qT>0.05)&&(qT<0.25)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackEta)<0.9)&&(invariantMassKmu<0.6)) {
360 if((kinkAngle>maxDecAngpimu)&&(qT>0.04)&&(qT<0.30)&&((kink->GetR()>120.)&&(kink->GetR()<220.))&&(TMath::Abs(trackEta)<0.9)&&(invariantMassKmu<0.6)) {
361
362         if( (kinkAngle>maxDecAngKmu*1.1) && ( motherMfromKink.Mag()> 1.2 ) ) fcodeH->Fill(TMath::Abs(code1), TMath::Abs(dcode1));
363                 if ( (kinkAngle>maxDecAngKmu*1.1) && ( motherMfromKink.Mag()> 1.2 ) )  fAngMomKC->Fill(motherMfromKink.Mag(), kinkAngle);
364
365
366                  if ( (kinkAngle>maxDecAngKmu*1.1) && ( motherMfromKink.Mag()> 1.2 ) ) continue; // maximum angle selection revised 22/10/2009
367
368                                     fHistPtKaon->Fill(track->Pt());   //all PID kink-kaon
369
370 //  kaons from k to mun and k to pipi and to e decay 
371          if(       ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==13))||    
372            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==11))  ||    
373            ( (TMath::Abs(code1)==321)&&(TMath::Abs(dcode1)==211))  ) { 
374
375          if ( label<=nPrim ) fPtPrKink->Fill(track->Pt());
376              fKinkKaon->Fill(track->Pt());        
377                 fHistEtaK->Fill(trackEta);
378                              }
379          else {
380              fKinkKaonBg->Fill(track->Pt());     
381  fdcodeH->Fill( TMath::Abs(code1), TMath::Abs(dcode1));   // put it here,  22/10/2009
382           }   // primary and all +BG    
383
384         }  //  kink selection 
385                   
386
387         }  //End Kink Information    
388   
389
390   } //track loop 
391
392   // loop over mc particles
393
394   // variable to count tracks
395   Int_t nMCTracks = 0;
396   Int_t nMCKinkKs = 0;
397
398 for (Int_t iMc = 0; iMc < nPrim; iMc++)
399   {
400    // AliDebug(AliLog::kDebug+1, Form("MC Loop: Processing particle %d.", iMc));
401
402     TParticle* particle = stack->Particle(iMc);
403
404     if (!particle)
405     {
406     //  AliDebug(AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack (mc loop).", iMc));
407       continue;
408     }
409
410  //   
411     Float_t eta = particle->Eta();
412
413     if (TMath::Abs(eta) > 0.9)
414       continue;
415
416             Double_t ptK = particle->Pt();
417
418    if( ptK <0.300) continue;
419
420       Int_t code = particle->GetPdgCode();
421             Int_t  mcProcess=-1011;
422       if ((code==321)||(code==-321)){
423             
424
425               fptKMC->Fill(ptK);
426
427     
428   // fMultiplMC->Fill(nPrim);
429        Int_t nMCKpi =0;
430
431             Int_t firstD=particle->GetFirstDaughter();
432             Int_t lastD=particle->GetLastDaughter();
433 //loop on secondaries
434              for (Int_t k=firstD;k<=lastD;k++) {
435               if ( k > 0 )    {
436              TParticle*    daughter1=stack->Particle(k);   // 27/8   
437              Int_t dcode = daughter1->GetPdgCode();
438
439         mcProcess=daughter1->GetUniqueID();
440      if (mcProcess==4) {        
441     
442
443                  // frad->Fill(daughter1->R());
444
445                  if (((daughter1->R())>120)&&((daughter1->R())<220) ){
446         if (( ( code==321 )&& ( dcode ==-13  ))|| (( code == -321 )&& ( dcode == 13))) fgenPtEtR->Fill( ptK );//to muon
447       //  if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211))) fgenPtEtR->Fill( ptK );//to pion
448         if (( ( code==321 )&& ( dcode ==-11  ))|| (( code == -321 )&& ( dcode ==11))) fgenPtEtR->Fill( ptK );// to electr
449         if (( ( code==321 )&& ( dcode ==211  ))|| (( code == -321 )&& ( dcode ==-211)))    nMCKpi++ ; 
450  
451   //fMultMCK->Fill(nPrim);
452                  frad->Fill(daughter1->R());
453     nMCKinkKs++;
454        }
455 //  next only  to mu decay
456                  if (((code==321)&&(dcode==-13))||((code==-321)&&(dcode==13))){
457
458                  
459                     if (((daughter1->R())>120)&&((daughter1->R())<220)){
460                  //   fgenpt->Fill(ptK);
461                     }
462                     }
463                     }//    decay
464          } // positive k
465        }//  daughters
466               if( nMCKpi == 1) fgenPtEtR->Fill(ptK);  //  k to pipi
467               if( nMCKpi > 1) fgenpt->Fill(ptK);// k to pipipi
468
469       }   /// kaons
470
471
472       
473     nMCTracks++;
474   }// end of mc particle
475
476 //  printf("hello mc");
477   fMultiplMC->Fill(nMCTracks);
478   if( nMCKinkKs>0) fMultMCK->Fill(nPrim);
479
480   PostData(1, fListOfHistos);
481
482 }      
483
484 //________________________________________________________________________
485 void AliAnalysisKinkESDMC::Terminate(Option_t *) 
486 {
487   // Draw result to the screen 
488   // Called once at the end of the query
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 }