PWG2/KINK/AliResonanceKink.cxx

   1 /**************************************************************************
   2  * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr)   *
   3  * Contributors are mentioned in the code where appropriate.              *
   4  *                                                                        *
   5  * Permission to use, copy, modify and distribute this software and its   *
   6  * documentation strictly for non-commercial purposes is hereby granted   *
   7  * without fee, provided that the above copyright notice appears in all   *
   8  * copies and that both the copyright notice and this permission notice   *
   9  * appear in the supporting documentation. The authors make no claims     *
  10  * about the suitability of this software for any purpose. It is          *
  11  * provided "as is" without express or implied warranty.                  *
  12  **************************************************************************/
  13
  14 //----------------------------------------------------------------------------------------------------------------
  15 //                        class AliResonanceKink
  16 //        Example of an analysis task for reconstructing resonances having at least one kaon-kink in their decay
  17 //        products. It provides basic plots as well as plots helping to calculate the corrections.
  18 //        Usage: To analyse a resonance having a kaon in its decay products, one has to modify the integer
  19 //        variables resonancePDG, fdaughter1pdg and fdaughter2pdg accordingly as well as daughter1pdgMass  and daughter2pdgMass.
  20 //        Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed
  21 //-----------------------------------------------------------------------------------------------------------------
  22
  23 #include "TChain.h"
  24 #include "TTree.h"
  25 #include "TH2D.h"
  26 #include "TParticle.h"
  27 #include <TVector3.h>
  28 #include <TDatabasePDG.h>
  29 #include <TParticlePDG.h>
  30 #include "TF1.h"
  31 #include "TList.h"
  32 #include "TString.h"
  33 #include "AliMCEventHandler.h"
  34 #include "AliMCEvent.h"
  35 #include "AliResonanceKink.h"
  36 #include "AliESDkink.h"
  37 #include "AliStack.h"
  38 #include "AliESDtrack.h"
  39 #include "AliESDEvent.h"
  40
  41 ClassImp(AliResonanceKink)
  42
  43 //________________________________________________________________________
  44 AliResonanceKink::AliResonanceKink()
  45   : TObject(), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(0),
  46   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(0), fLowX(0), fHighX(0), fdaughter1pdg(0), fdaughter2pdg(0), fresonancePDGcode(0)
  47
  48 {
  49   // Constructor
  50 }
  51
  52 //________________________________________________________________________
  53 AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG)
  54   : TObject(), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(0),
  55   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG)
  56
  57 {
  58    // Constructor
  59
  60    fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
  61
  62    fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
  63    fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
  64    fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
  65
  66    fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
  67    fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
  68    fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
  69    fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
  70    fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
  71    fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
  72    fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
  73    fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
  74    fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
  75    fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
  76    fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
  77
  78    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);
  79    f1->SetParameter(0,0.493677);
  80    f1->SetParameter(1,0.9127037);
  81    f1->SetParameter(2,TMath::Pi());
  82
  83    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);
  84    f2->SetParameter(0,0.13957018);
  85    f2->SetParameter(1,0.2731374);
  86    f2->SetParameter(2,TMath::Pi());
  87
  88    fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
  89
  90 }
  91
  92 //________________________________________________________________________
  93 AliResonanceKink:: ~AliResonanceKink()
  94 {
  95  //  Destructor
  96  if(fOpeningAngle) delete fOpeningAngle;
  97  if(fInvariantMass) delete fInvariantMass;
  98  if(fInvMassTrue) delete fInvMassTrue;
  99  if(fPhiBothKinks) delete fPhiBothKinks;
 100  if(fRecPt) delete fRecPt;
 101  if(fRecEta) delete fRecEta;
 102  if(fRecEtaPt) delete fRecEtaPt;
 103  if(fSimPt) delete fSimPt;
 104  if(fSimEta) delete fSimEta;
 105  if(fSimEtaPt) delete fSimEtaPt;
 106  if(fSimPtKink) delete fSimPtKink;
 107  if(fSimEtaKink) delete fSimEtaKink;
 108  if(fSimEtaPtKink) delete fSimEtaPtKink;
 109  if(fhdr) delete fhdr;
 110  if(fhdz) delete fhdz;
 111  if(fvtxz) delete fvtxz;
 112 }
 113 //________________________________________________________________________
 114 TList* AliResonanceKink::GetHistogramList()
 115 {
 116   // Adding histograms to the list
 117   fListOfHistos=new TList();
 118
 119   fListOfHistos->Add(fOpeningAngle);
 120   fListOfHistos->Add(fInvariantMass);
 121   fListOfHistos->Add(fInvMassTrue);
 122   fListOfHistos->Add(fPhiBothKinks);
 123   fListOfHistos->Add(fRecPt);
 124   fListOfHistos->Add(fRecEta);
 125   fListOfHistos->Add(fRecEtaPt);
 126   fListOfHistos->Add(fSimPt);
 127   fListOfHistos->Add(fSimEta);
 128   fListOfHistos->Add(fSimEtaPt);
 129   fListOfHistos->Add(fSimPtKink);
 130   fListOfHistos->Add(fSimEtaKink);
 131   fListOfHistos->Add(fSimEtaPtKink);
 132   fListOfHistos->Add(fhdr);
 133   fListOfHistos->Add(fhdz);
 134   fListOfHistos->Add(fvtxz);
 135
 136   return fListOfHistos;
 137 }
 138
 139 //________________________________________________________________________
 140 void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
 141 {
 142   //  Initialisation of the output histograms
 143   fNbins=nbins;
 144   fLowX=nlowx;
 145   fHighX=nhighx;
 146
 147   fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
 148
 149   fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
 150   fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
 151   fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
 152
 153   fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
 154   fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
 155   fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 156   fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
 157   fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
 158   fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 159   fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
 160   fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
 161   fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
 162   fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
 163   fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
 164
 165   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);
 166   f1->SetParameter(0,0.493677);
 167   f1->SetParameter(1,0.9127037);
 168   f1->SetParameter(2,TMath::Pi());
 169
 170   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);
 171   f2->SetParameter(0,0.13957018);
 172   f2->SetParameter(1,0.2731374);
 173   f2->SetParameter(2,TMath::Pi());
 174
 175   fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
 176 }
 177
 178 //________________________________________________________________________
 179 void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
 180 {
 181   // Main loop
 182   // Called for each event
 183   Int_t resonancePDGcode, antiresonancePDGcode;
 184
 185   if (fdaughter1pdg==kdaughterKaon)  {
 186     resonancePDGcode=fresonancePDGcode;
 187     antiresonancePDGcode=-fresonancePDGcode;
 188   }
 189   if (fdaughter1pdg!=kdaughterKaon)  {
 190     resonancePDGcode=-fresonancePDGcode;
 191     antiresonancePDGcode=fresonancePDGcode;
 192   }
 193   if (fdaughter1pdg==fdaughter2pdg)  {
 194     resonancePDGcode=fresonancePDGcode;
 195     antiresonancePDGcode=fresonancePDGcode;
 196   }
 197
 198    Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
 199    Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
 200
 201   if (!esd) {
 202     Printf("ERROR: fESD not available");
 203     return;
 204   }
 205
 206     if (!mcEvent) {
 207     Printf("ERROR: mcEvent not available");
 208     return;
 209   }
 210
 211   AliStack* stack=mcEvent->Stack();
 212
 213   if(fAnalysisType == "MC") {
 214   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
 215   {
 216     TParticle* particle = stack->Particle(iMc);
 217
 218     if (!particle)
 219     {
 220       Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
 221       continue;
 222     }
 223
 224      if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
 225        Int_t firstD=particle->GetFirstDaughter();
 226        Int_t lastD=particle->GetLastDaughter();
 227        TParticle *daughterParticle1=stack->Particle(firstD);
 228        TParticle *daughterParticle2=stack->Particle(lastD);
 229
 230        TParticle* kaonFirstDaughter;
 231        Int_t mcProcessKaonFirstDaughter = -999;
 232
 233        for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++){
 234         if ((daughterParticle1->GetFirstDaughter()+ia)!=-1) {
 235           kaonFirstDaughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
 236           mcProcessKaonFirstDaughter=kaonFirstDaughter->GetUniqueID();
 237         }
 238        }
 239
 240        if((daughterParticle1->Pt()>0.25)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<1.1)&&            (TMath::Abs(daughterParticle2->Eta())<1.1)&&(TMath::Abs(particle->Eta())<1.1)) {
 241          fSimEta->Fill(particle->Eta());
 242          fSimPt->Fill(particle->Pt());
 243          fSimEtaPt->Fill(particle->Pt(), particle->Eta());
 244          if(mcProcessKaonFirstDaughter==4) {
 245            fSimPtKink->Fill(particle->Pt());
 246            fSimEtaKink->Fill(particle->Eta());
 247            fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
 248          }
 249        }
 250      }
 251   }
 252
 253   } // end fAnalysisType==MC
 254   else
 255
 256   if(fAnalysisType == "ESD") {
 257   const AliESDVertex* vertex = GetEventVertex(esd);
 258   if(!vertex) return;
 259   Double_t vtx[3];
 260   vertex->GetXYZ(vtx);
 261   fvtxz->Fill(vtx[2]);
 262   Double_t ptrackpos[3], ptrackneg[3];
 263
 264   TLorentzVector p4pos, anp4pos;
 265   TLorentzVector p4neg, anp4neg;
 266   TLorentzVector p4comb, anp4comb;
 267
 268   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 269     AliESDtrack* trackpos = esd->GetTrack(iTracks);
 270     if (!trackpos) {
 271       Printf("ERROR: Could not receive track %d", iTracks);
 272       continue;
 273     }
 274     if (trackpos->GetSign() < 0) continue;
 275
 276     trackpos->GetPxPyPz(ptrackpos);
 277
 278     Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
 279
 280     Float_t bpos[2];
 281     Float_t bCovpos[3];
 282     trackpos->GetImpactParameters(bpos,bCovpos);
 283
 284     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 285      Printf("Estimated b resolution lower or equal zero!");
 286      bCovpos[0]=0; bCovpos[2]=0;
 287     }
 288
 289     Float_t dcaToVertexXYpos = bpos[0];
 290     Float_t dcaToVertexZpos = bpos[1];
 291
 292     fhdr->Fill(dcaToVertexXYpos);
 293     fhdz->Fill(dcaToVertexZpos);
 294
 295     if(nSigmaToVertex>=4) continue;
 296     if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
 297
 298     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
 299
 300     if(posTrackMom.Perp()<=0.25) continue;
 301
 302     TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
 303     if (!partpos) continue;
 304     Int_t pdgpos = partpos->GetPdgCode();
 305     Int_t mumlabelpos=partpos->GetFirstMother();
 306     mumlabelpos = TMath::Abs(mumlabelpos);
 307     TParticle * mumpos=stack->Particle(mumlabelpos);
 308     if (!mumpos) continue;
 309     Int_t mumpdgpos = mumpos->GetPdgCode();
 310
 311     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
 312     Int_t kaonKinkFlag=0;
 313     if(indexKinkPos<0){
 314
 315         AliESDkink *poskink=esd->GetKink(TMath::Abs(indexKinkPos)-1);
 316         const TVector3 motherMfromKinkPos(poskink->GetMotherP());
 317         const TVector3 daughterMKinkPos(poskink->GetDaughterP());
 318         Float_t posQt=poskink->GetQt();
 319
 320         Double_t maxDecAngKmuPos=f1->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
 321         Double_t maxDecAngpimuPos=f2->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
 322
 323         Float_t kinkAnglePos=TMath::RadToDeg()*poskink->GetAngle(2);
 324
 325         Float_t energyDaughterMu=TMath::Sqrt(daughterMKinkPos.Mag()*daughterMKinkPos.Mag()+0.105658*0.105658);
 326         Float_t p1XM= motherMfromKinkPos.Px();
 327         Float_t p1YM= motherMfromKinkPos.Py();
 328         Float_t p1ZM= motherMfromKinkPos.Pz();
 329         Float_t p2XM= daughterMKinkPos.Px();
 330         Float_t p2YM= daughterMKinkPos.Py();
 331         Float_t p2ZM= daughterMKinkPos.Pz();
 332         Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 333         Double_t invariantMassKmuPos= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKinkPos.Mag()*motherMfromKinkPos.Mag());
 334
 335         if((kinkAnglePos>maxDecAngpimuPos)&&(posQt>0.05)&&(posQt<0.25)&&((poskink->GetR()>110.)&&(poskink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmuPos<0.6)) {
 336
 337           if(posTrackMom.Mag()<=1.1) {
 338            kaonKinkFlag=1;
 339           }
 340           else
 341           if (kinkAnglePos<maxDecAngKmuPos) {
 342            kaonKinkFlag=1;
 343           }
 344         }
 345
 346     }  //End Kink Information
 347
 348     if(kaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
 349
 350     if(indexKinkPos==0) {
 351       UInt_t status=trackpos->GetStatus();
 352       if((status&AliESDtrack::kTPCrefit)==0) continue;
 353       if(trackpos->GetTPCclusters(0)<50) continue;
 354       if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
 355       Double_t extCovPos[15];
 356       trackpos->GetExternalCovariance(extCovPos);
 357       if(extCovPos[0]>2) continue;
 358       if(extCovPos[2]>2) continue;
 359       if(extCovPos[5]>0.5) continue;
 360       if(extCovPos[9]>0.5) continue;
 361       if(extCovPos[14]>2) continue;
 362
 363       p4pos.SetVectM(posTrackMom, daughter2pdgMass);
 364
 365     }
 366
 367       for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
 368         if(iTracks==j) continue;
 369         AliESDtrack* trackneg=esd->GetTrack(j);
 370         if (trackneg->GetSign() > 0) continue;
 371
 372         trackneg->GetPxPyPz(ptrackneg);
 373         Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
 374
 375         Float_t bneg[2];
 376         Float_t bCovneg[3];
 377         trackneg->GetImpactParameters(bneg,bCovneg);
 378         if (bCovneg[0]<=0 || bCovneg[2]<=0) {
 379           Printf("Estimated b resolution lower or equal zero!");
 380           bCovneg[0]=0; bCovneg[2]=0;
 381         }
 382
 383         Float_t dcaToVertexXYneg = bneg[0];
 384         Float_t dcaToVertexZneg = bneg[1];
 385
 386         fhdr->Fill(dcaToVertexXYneg);
 387         fhdz->Fill(dcaToVertexZneg);
 388
 389         if(negSigmaToVertex>=4) continue;
 390         if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
 391
 392         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
 393
 394         if(negTrackMom.Perp()<=0.25) continue;
 395
 396         TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
 397         if (!partneg) continue;
 398         Int_t pdgneg = partneg->GetPdgCode();
 399         Int_t mumlabelneg=partneg->GetFirstMother();
 400         mumlabelneg = TMath::Abs(mumlabelneg);
 401         TParticle * mumneg=stack->Particle(mumlabelneg);
 402         if (!mumneg) continue;
 403         Int_t mumpdgneg = mumneg->GetPdgCode();
 404
 405         Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
 406         Int_t negKaonKinkFlag=0;
 407         if(indexKinkNeg<0){
 408
 409          AliESDkink *negkink=esd->GetKink(TMath::Abs(indexKinkNeg)-1);
 410          const TVector3 motherMfromKinkNeg(negkink->GetMotherP());
 411          const TVector3 daughterMKinkNeg(negkink->GetDaughterP());
 412          Float_t negQt=negkink->GetQt();
 413
 414          Double_t maxDecAngKmuNeg=f1->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
 415          Double_t maxDecAngpimuNeg=f2->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
 416
 417          Float_t kinkAngleNeg=TMath::RadToDeg()*negkink->GetAngle(2);
 418
 419          Float_t energyDaughterMuNeg=TMath::Sqrt(daughterMKinkNeg.Mag()*daughterMKinkNeg.Mag()+0.105658*0.105658);
 420          Float_t p1XMNeg= motherMfromKinkNeg.Px();
 421          Float_t p1YMNeg= motherMfromKinkNeg.Py();
 422          Float_t p1ZMNeg= motherMfromKinkNeg.Pz();
 423          Float_t p2XMNeg= daughterMKinkNeg.Px();
 424          Float_t p2YMNeg= daughterMKinkNeg.Py();
 425          Float_t p2ZMNeg= daughterMKinkNeg.Pz();
 426          Float_t p3DaughterNeg=TMath::Sqrt(((p1XMNeg-p2XMNeg)*(p1XMNeg-p2XMNeg))+((p1YMNeg-p2YMNeg)*(p1YMNeg-p2YMNeg))+((p1ZMNeg-p2ZMNeg)*(p1ZMNeg-p2ZMNeg)));
 427          Double_t invariantMassKmuNeg= TMath::Sqrt((energyDaughterMuNeg+p3DaughterNeg)*(energyDaughterMuNeg+p3DaughterNeg)-motherMfromKinkNeg.Mag()*motherMfromKinkNeg.Mag());
 428
 429          if((kinkAngleNeg>maxDecAngpimuNeg)&&(negQt>0.05)&&(negQt<0.25)&&((negkink->GetR()>110.)&&(negkink->GetR()<230.))&&(TMath::Abs(negTrackMom.Eta())<1.1)&&(invariantMassKmuNeg<0.6)) {
 430
 431            if(negTrackMom.Mag()<=1.1) {
 432                 negKaonKinkFlag=1;
 433            }
 434            else
 435            if (kinkAngleNeg<maxDecAngKmuNeg) {
 436                 negKaonKinkFlag=1;
 437            }
 438          }
 439
 440         }  //End Kink Information
 441
 442         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
 443
 444         if(indexKinkNeg==0)  {
 445            UInt_t statusneg=trackneg->GetStatus();
 446
 447            if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
 448
 449            if(trackneg->GetTPCclusters(0)<50) continue;
 450            if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
 451            Double_t extCovneg[15];
 452            trackneg->GetExternalCovariance(extCovneg);
 453            if(extCovneg[0]>2) continue;
 454            if(extCovneg[2]>2) continue;
 455            if(extCovneg[5]>0.5) continue;
 456            if(extCovneg[9]>0.5) continue;
 457            if(extCovneg[14]>2) continue;
 458
 459           anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
 460
 461         }
 462
 463         Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
 464
 465         if((kaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
 466          p4comb=anp4pos;
 467          p4comb+=p4neg;
 468          if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
 469         }
 470
 471         if(negKaonKinkFlag==1) {
 472           p4comb=p4pos;
 473           p4comb+=p4neg;
 474           fInvariantMass->Fill(p4comb.M());
 475           if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
 476           &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
 477             fOpeningAngle->Fill(openingAngle);
 478             fInvMassTrue->Fill(p4comb.M());
 479             if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
 480               fRecPt->Fill(p4comb.Vect().Pt());
 481               fRecEta->Fill(p4comb.Vect().Eta());
 482               fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
 483
 484             }
 485
 486            }
 487
 488         }
 489
 490         if(kaonKinkFlag==1) {
 491           anp4comb=anp4pos;
 492           anp4comb+=anp4neg;
 493           fInvariantMass->Fill(anp4comb.M());
 494           if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
 495           &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
 496             fOpeningAngle->Fill(openingAngle);
 497             fInvMassTrue->Fill(p4comb.M());
 498             if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
 499              fRecPt->Fill(anp4comb.Vect().Pt());
 500              fRecEta->Fill(anp4comb.Vect().Eta());
 501              fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
 502            }
 503
 504          }
 505
 506         }
 507
 508       } //inner track loop
 509
 510   } //outer track loop
 511
 512   } // end fAnalysisType == ESD
 513
 514 }
 515
 516 //____________________________________________________________________//
 517 Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 518   // Calculates the number of sigma to the vertex.
 519
 520   Float_t b[2];
 521   Float_t bRes[2];
 522   Float_t bCov[3];
 523
 524     esdTrack->GetImpactParameters(b,bCov);
 525
 526   if (bCov[0]<=0 || bCov[2]<=0) {
 527     bCov[0]=0; bCov[2]=0;
 528   }
 529
 530   bRes[0] = TMath::Sqrt(bCov[0]);
 531   bRes[1] = TMath::Sqrt(bCov[2]);
 532
 533   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 534
 535   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 536
 537   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 538
 539   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 540
 541   return d;
 542 }
 543
 544 //________________________________________________________________________
 545 const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
 546 {
 547   // Get the vertex
 548
 549   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 550
 551   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 552   else
 553   {
 554      vertex = esd->GetPrimaryVertexSPD();
 555      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 556      else
 557      return 0;
 558   }
 559 }
 560
 561
 562