PWG2/KINK/AliResonanceKinkPID.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 AliResonanceKinkPID
  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, daughter1 and daughter2 accordingly as well as daughter1Mass  and daughter2Mass.
  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 "AliAnalysisTaskSE.h"
  32 #include "AliAnalysisManager.h"
  33
  34 #include "AliESDInputHandler.h"
  35
  36 #include "AliMCEventHandler.h"
  37 #include "AliMCEvent.h"
  38
  39 #include "AliResonanceKinkPID.h"
  40 #include "AliESDkink.h"
  41 #include "AliStack.h"
  42
  43 ClassImp(AliResonanceKinkPID)
  44
  45 //________________________________________________________________________
  46 AliResonanceKinkPID::AliResonanceKinkPID()
  47   : AliAnalysisTaskSE(), fESD(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0),
  48   fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(0),
  49   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType("ESD"), fvtxz(0)
  50
  51 {
  52   // Constructor
  53 }
  54
  55 //________________________________________________________________________
  56 AliResonanceKinkPID::AliResonanceKinkPID(const char *name)
  57   : AliAnalysisTaskSE(name), fESD(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fRecPt(0), fRecEta(0), fRecEtaPt(0),
  58   fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(0),
  59   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType("ESD"), fvtxz(0)
  60
  61 {
  62   // Constructor
  63
  64   // Define input and output slots here
  65   // Input slot #0 works with a TChain
  66   DefineInput(0, TChain::Class());
  67   DefineOutput(1, TList::Class());
  68 }
  69
  70 //________________________________________________________________________
  71 void AliResonanceKinkPID::ConnectInputData(Option_t *)
  72 {
  73   // Connect ESD or AOD here
  74   // Called once
  75
  76   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  77   if (!tree) {
  78     Printf("ERROR: Could not read chain from input slot 0");
  79   } else {
  80
  81     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  82
  83     if (!esdH) {
  84       Printf("ERROR: Could not get ESDInputHandler");
  85     } else
  86       fESD = esdH->GetEvent();
  87   }
  88 }
  89
  90 //________________________________________________________________________
  91 void AliResonanceKinkPID::CreateOutputObjects()
  92 {
  93   // Create histograms
  94   // Called once
  95
  96    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);
  97    f1->SetParameter(0,0.493677);
  98    f1->SetParameter(1,0.9127037);
  99    f1->SetParameter(2,TMath::Pi());
 100
 101    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);
 102    f2->SetParameter(0,0.13957018);
 103    f2->SetParameter(1,0.2731374);
 104    f2->SetParameter(2,TMath::Pi());
 105
 106    // OpenFile(1); Uncomment for proof analysis
 107    fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
 108
 109    // for K*0(892)
 110    fInvariantMass=new TH1D("fInvariantMass"," ", 60,0.6,1.2);
 111    fInvMassTrue=new TH1D("fInvMassTrue"," ",60,0.6,1.2);
 112    fPhiBothKinks=new TH1D("fPhiBothKinks"," ", 60,0.6,1.2);  // Not applicable for K*0
 113
 114    //for phi(1020)
 115    //fInvariantMass=new TH1D("fInvariantMass"," ", 70,0.99,1.088);
 116    //fInvMassTrue=new TH1D("fInvMassTrue"," ", 70,0.99,1.088);
 117    //fPhiBothKinks=new TH1D("fPhiBothKinks"," ", 70,0.99,1.088);
 118
 119    fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
 120    fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
 121    fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 122    fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
 123    fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
 124    fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 125    fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
 126    fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
 127    fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
 128    fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
 129    fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
 130    fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
 131
 132    fListOfHistos=new TList();
 133
 134    fListOfHistos->Add(fOpeningAngle);
 135    fListOfHistos->Add(fInvariantMass);
 136    fListOfHistos->Add(fInvMassTrue);
 137    fListOfHistos->Add(fPhiBothKinks);
 138    fListOfHistos->Add(fRecPt);
 139    fListOfHistos->Add(fRecEta);
 140    fListOfHistos->Add(fRecEtaPt);
 141    fListOfHistos->Add(fSimPt);
 142    fListOfHistos->Add(fSimEta);
 143    fListOfHistos->Add(fSimEtaPt);
 144    fListOfHistos->Add(fSimPtKink);
 145    fListOfHistos->Add(fSimEtaKink);
 146    fListOfHistos->Add(fSimEtaPtKink);
 147    fListOfHistos->Add(fhdr);
 148    fListOfHistos->Add(fhdz);
 149    fListOfHistos->Add(fvtxz);
 150
 151 }
 152
 153 //________________________________________________________________________
 154 void AliResonanceKinkPID::Exec(Option_t *)
 155 {
 156   // Main loop
 157   // Called for each event
 158
 159    enum ResonanceType {kPhi=333, kKstar0=313, kLambda1520=3124};
 160    enum DaughterType {kdaughterPion=211, kdaughterKaon=321, kdaughterProton=2212};
 161
 162    Int_t resonancePDG=kKstar0;
 163
 164    Int_t daughter1=kdaughterKaon;
 165    Int_t daughter2=kdaughterPion;
 166
 167    Double_t daughter1Mass=TDatabasePDG::Instance()->GetParticle(kKPlus)->Mass();
 168    Double_t daughter2Mass=TDatabasePDG::Instance()->GetParticle(kPiPlus)->Mass();
 169
 170    AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 171    if (!eventHandler) {
 172     Printf("ERROR: Could not retrieve MC event handler");
 173     return;
 174    }
 175
 176   AliMCEvent* mcEvent = eventHandler->MCEvent();
 177   if (!mcEvent) {
 178      Printf("ERROR: Could not retrieve MC event");
 179      return;
 180   }
 181
 182   Printf("MC particles: %d", mcEvent->GetNumberOfTracks());
 183
 184   if (!fESD) {
 185     Printf("ERROR: fESD not available");
 186     return;
 187   }
 188
 189   AliStack* stack=mcEvent->Stack();
 190
 191   if(fAnalysisType == "MC") {
 192
 193   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
 194   {
 195     TParticle* particle = stack->Particle(iMc);
 196
 197     if (!particle)
 198     {
 199       Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
 200       continue;
 201     }
 202
 203      if(TMath::Abs(particle->GetPdgCode())==resonancePDG) {
 204        Int_t firstD=particle->GetFirstDaughter();
 205        Int_t lastD=particle->GetLastDaughter();
 206        TParticle *daughterParticle1=stack->Particle(firstD);
 207        TParticle *daughterParticle2=stack->Particle(lastD);
 208
 209        TParticle* kaonFirstDaughter;
 210        Int_t mcProcessKaonFirstDaughter;
 211
 212        for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++){
 213         if ((daughterParticle1->GetFirstDaughter()+ia)!=-1) {
 214           kaonFirstDaughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
 215           mcProcessKaonFirstDaughter=kaonFirstDaughter->GetUniqueID();
 216         }
 217        }
 218
 219        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)) {
 220          fSimEta->Fill(particle->Eta());
 221          fSimPt->Fill(particle->Pt());
 222          fSimEtaPt->Fill(particle->Pt(), particle->Eta());
 223          if(mcProcessKaonFirstDaughter==4) {
 224            fSimPtKink->Fill(particle->Pt());
 225            fSimEtaKink->Fill(particle->Eta());
 226            fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
 227          }
 228        }
 229      }
 230   }
 231
 232   } // end fAnalysisType==MC
 233   else
 234
 235   if(fAnalysisType == "ESD") {
 236
 237   const AliESDVertex* vertex = GetEventVertex(fESD);
 238   if(!vertex) return;
 239   Double_t vtx[3];
 240   vertex->GetXYZ(vtx);
 241   fvtxz->Fill(vtx[2]);
 242
 243   Double_t ptrackpos[3], ptrackneg[3];
 244
 245   TLorentzVector p4pos, anp4pos;
 246   TLorentzVector p4neg, anp4neg;
 247   TLorentzVector p4comb, anp4comb;
 248
 249
 250   for (Int_t iTracks = 0; iTracks < fESD->GetNumberOfTracks(); iTracks++) {
 251     AliESDtrack* trackpos = fESD->GetTrack(iTracks);
 252     if (!trackpos) {
 253       Printf("ERROR: Could not receive track %d", iTracks);
 254       continue;
 255     }
 256     if (trackpos->GetSign() < 0) continue;
 257
 258     trackpos->GetPxPyPz(ptrackpos);
 259
 260     Float_t nSigmaToVertex = GetSigmaToVertex(trackpos);
 261
 262     Float_t bpos[2];
 263     Float_t bCovpos[3];
 264     trackpos->GetImpactParameters(bpos,bCovpos);
 265
 266     if (bCovpos[0]<=0 || bCovpos[2]<=0) {
 267      Printf("Estimated b resolution lower or equal zero!");
 268      bCovpos[0]=0; bCovpos[2]=0;
 269     }
 270
 271     Float_t dcaToVertexXYpos = bpos[0];
 272     Float_t dcaToVertexZpos = bpos[1];
 273
 274     fhdr->Fill(dcaToVertexXYpos);
 275     fhdz->Fill(dcaToVertexZpos);
 276
 277     if(nSigmaToVertex>=4) continue;
 278     if((dcaToVertexXYpos>3.0)||(dcaToVertexZpos>3.0)) continue;
 279
 280     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
 281
 282     if(posTrackMom.Perp()<=0.25) continue;
 283
 284     TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
 285     if (!partpos) continue;
 286     Int_t pdgpos = partpos->GetPdgCode();
 287     Int_t mumlabelpos=partpos->GetFirstMother();
 288     mumlabelpos = TMath::Abs(mumlabelpos);
 289     TParticle * mumpos=stack->Particle(mumlabelpos);
 290     if (!mumpos) continue;
 291     Int_t mumpdgpos = mumpos->GetPdgCode();
 292
 293     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
 294     Int_t kaonKinkFlag=0;
 295     if(indexKinkPos<0){
 296
 297         AliESDkink *poskink=fESD->GetKink(TMath::Abs(indexKinkPos)-1);
 298         const TVector3 motherMfromKinkPos(poskink->GetMotherP());
 299         const TVector3 daughterMKinkPos(poskink->GetDaughterP());
 300         Float_t posQt=poskink->GetQt();
 301
 302         Double_t maxDecAngKmuPos=f1->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
 303         Double_t maxDecAngpimuPos=f2->Eval(motherMfromKinkPos.Mag(),0.,0.,0.);
 304
 305         Float_t kinkAnglePos=TMath::RadToDeg()*poskink->GetAngle(2);
 306
 307         Float_t energyDaughterMu=TMath::Sqrt(daughterMKinkPos.Mag()*daughterMKinkPos.Mag()+0.105658*0.105658);
 308         Float_t p1XM= motherMfromKinkPos.Px();
 309         Float_t p1YM= motherMfromKinkPos.Py();
 310         Float_t p1ZM= motherMfromKinkPos.Pz();
 311         Float_t p2XM= daughterMKinkPos.Px();
 312         Float_t p2YM= daughterMKinkPos.Py();
 313         Float_t p2ZM= daughterMKinkPos.Pz();
 314         Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 315         Double_t invariantMassKmuPos= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKinkPos.Mag()*motherMfromKinkPos.Mag());
 316
 317         if((kinkAnglePos>maxDecAngpimuPos)&&(posQt>0.05)&&(posQt<0.25)&&((poskink->GetR()>110.)&&(poskink->GetR()<230.))&&(TMath::Abs(posTrackMom.Eta())<1.1)&&(invariantMassKmuPos<0.6)) {
 318
 319           if(posTrackMom.Mag()<=1.1) {
 320            kaonKinkFlag=1;
 321           }
 322           else
 323           if (kinkAnglePos<maxDecAngKmuPos) {
 324            kaonKinkFlag=1;
 325           }
 326         }
 327
 328     }  //End Kink Information
 329
 330     if(kaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1Mass);
 331
 332     if(indexKinkPos==0) {
 333       UInt_t status=trackpos->GetStatus();
 334       if((status&AliESDtrack::kTPCrefit)==0) continue;
 335       if(trackpos->GetTPCclusters(0)<50) continue;
 336       if((trackpos->GetTPCchi2()/trackpos->GetTPCclusters(0))>3.5) continue;
 337       Double_t extCovPos[15];
 338       trackpos->GetExternalCovariance(extCovPos);
 339       if(extCovPos[0]>2) continue;
 340       if(extCovPos[2]>2) continue;
 341       if(extCovPos[5]>0.5) continue;
 342       if(extCovPos[9]>0.5) continue;
 343       if(extCovPos[14]>2) continue;
 344
 345       p4pos.SetVectM(posTrackMom, daughter2Mass);
 346
 347     }
 348
 349       for (Int_t j=0; j<fESD->GetNumberOfTracks(); j++) {
 350         if(iTracks==j) continue;
 351         AliESDtrack* trackneg=fESD->GetTrack(j);
 352         if (trackneg->GetSign() > 0) continue;
 353
 354         trackneg->GetPxPyPz(ptrackneg);
 355         Float_t negSigmaToVertex = GetSigmaToVertex(trackneg);
 356
 357         Float_t bneg[2];
 358         Float_t bCovneg[3];
 359         trackneg->GetImpactParameters(bneg,bCovneg);
 360         if (bCovneg[0]<=0 || bCovneg[2]<=0) {
 361           Printf("Estimated b resolution lower or equal zero!");
 362           bCovneg[0]=0; bCovneg[2]=0;
 363         }
 364
 365         Float_t dcaToVertexXYneg = bneg[0];
 366         Float_t dcaToVertexZneg = bneg[1];
 367
 368         fhdr->Fill(dcaToVertexXYneg);
 369         fhdz->Fill(dcaToVertexZneg);
 370
 371         if(negSigmaToVertex>=4) continue;
 372         if((dcaToVertexXYneg>3.0)||(dcaToVertexZneg>3.0)) continue;
 373
 374         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
 375
 376         if(negTrackMom.Perp()<=0.25) continue;
 377
 378         TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
 379         if (!partneg) continue;
 380         Int_t pdgneg = partneg->GetPdgCode();
 381         Int_t mumlabelneg=partneg->GetFirstMother();
 382         mumlabelneg = TMath::Abs(mumlabelneg);
 383         TParticle * mumneg=stack->Particle(mumlabelneg);
 384         if (!mumneg) continue;
 385         Int_t mumpdgneg = mumneg->GetPdgCode();
 386
 387         Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
 388         Int_t negKaonKinkFlag=0;
 389         if(indexKinkNeg<0){
 390
 391          AliESDkink *negkink=fESD->GetKink(TMath::Abs(indexKinkNeg)-1);
 392          const TVector3 motherMfromKinkNeg(negkink->GetMotherP());
 393          const TVector3 daughterMKinkNeg(negkink->GetDaughterP());
 394          Float_t negQt=negkink->GetQt();
 395
 396          Double_t maxDecAngKmuNeg=f1->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
 397          Double_t maxDecAngpimuNeg=f2->Eval(motherMfromKinkNeg.Mag(),0.,0.,0.);
 398
 399          Float_t kinkAngleNeg=TMath::RadToDeg()*negkink->GetAngle(2);
 400
 401          Float_t energyDaughterMuNeg=TMath::Sqrt(daughterMKinkNeg.Mag()*daughterMKinkNeg.Mag()+0.105658*0.105658);
 402          Float_t p1XMNeg= motherMfromKinkNeg.Px();
 403          Float_t p1YMNeg= motherMfromKinkNeg.Py();
 404          Float_t p1ZMNeg= motherMfromKinkNeg.Pz();
 405          Float_t p2XMNeg= daughterMKinkNeg.Px();
 406          Float_t p2YMNeg= daughterMKinkNeg.Py();
 407          Float_t p2ZMNeg= daughterMKinkNeg.Pz();
 408          Float_t p3DaughterNeg=TMath::Sqrt(((p1XMNeg-p2XMNeg)*(p1XMNeg-p2XMNeg))+((p1YMNeg-p2YMNeg)*(p1YMNeg-p2YMNeg))+((p1ZMNeg-p2ZMNeg)*(p1ZMNeg-p2ZMNeg)));
 409          Double_t invariantMassKmuNeg= TMath::Sqrt((energyDaughterMuNeg+p3DaughterNeg)*(energyDaughterMuNeg+p3DaughterNeg)-motherMfromKinkNeg.Mag()*motherMfromKinkNeg.Mag());
 410
 411          if((kinkAngleNeg>maxDecAngpimuNeg)&&(negQt>0.05)&&(negQt<0.25)&&((negkink->GetR()>110.)&&(negkink->GetR()<230.))&&(TMath::Abs(negTrackMom.Eta())<1.1)&&(invariantMassKmuNeg<0.6)) {
 412
 413            if(negTrackMom.Mag()<=1.1) {
 414                 negKaonKinkFlag=1;
 415            }
 416            else
 417            if (kinkAngleNeg<maxDecAngKmuNeg) {
 418                 negKaonKinkFlag=1;
 419            }
 420          }
 421
 422         }  //End Kink Information
 423
 424         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1Mass);
 425
 426         if(indexKinkNeg==0)  {
 427            UInt_t statusneg=trackneg->GetStatus();
 428
 429            if((statusneg&AliESDtrack::kTPCrefit)==0) continue;
 430
 431            if(trackneg->GetTPCclusters(0)<50) continue;
 432            if((trackneg->GetTPCchi2()/trackneg->GetTPCclusters(0))>3.5) continue;
 433            Double_t extCovneg[15];
 434            trackneg->GetExternalCovariance(extCovneg);
 435            if(extCovneg[0]>2) continue;
 436            if(extCovneg[2]>2) continue;
 437            if(extCovneg[5]>0.5) continue;
 438            if(extCovneg[9]>0.5) continue;
 439            if(extCovneg[14]>2) continue;
 440
 441           anp4neg.SetVectM(negTrackMom, daughter2Mass);
 442
 443         }
 444
 445         Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
 446
 447         if((kaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
 448          p4comb=anp4pos;
 449          p4comb+=p4neg;
 450          if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
 451         }
 452
 453         if(negKaonKinkFlag==1) {
 454           p4comb=p4pos;
 455           p4comb+=p4neg;
 456           fInvariantMass->Fill(p4comb.M());
 457           if ((mumpdgpos==(-resonancePDG))&&(mumpdgneg==(-resonancePDG))&&(mumlabelpos==mumlabelneg)
 458           &&(pdgpos==daughter2)&&(pdgneg==(-daughter1))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
 459             fOpeningAngle->Fill(openingAngle);
 460             fInvMassTrue->Fill(p4comb.M());
 461             if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
 462               fRecPt->Fill(p4comb.Vect().Pt());
 463               fRecEta->Fill(p4comb.Vect().Eta());
 464               fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
 465
 466             }
 467
 468            }
 469
 470         }
 471
 472         if(kaonKinkFlag==1) {
 473           anp4comb=anp4pos;
 474           anp4comb+=anp4neg;
 475           fInvariantMass->Fill(anp4comb.M());
 476           if ((mumpdgpos==resonancePDG)&&(mumpdgneg==resonancePDG)&&(mumlabelpos==mumlabelneg)
 477           &&(pdgpos==daughter1)&&(pdgneg==(-daughter2))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
 478             fOpeningAngle->Fill(openingAngle);
 479             fInvMassTrue->Fill(p4comb.M());
 480             if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
 481              fRecPt->Fill(anp4comb.Vect().Pt());
 482              fRecEta->Fill(anp4comb.Vect().Eta());
 483              fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
 484            }
 485
 486          }
 487
 488         }
 489
 490       } //inner track loop
 491
 492   } //outer track loop
 493
 494   } // end fAnalysisType == ESD
 495
 496     PostData(1, fListOfHistos);
 497 }
 498
 499 //________________________________________________________________________
 500 void AliResonanceKinkPID::Terminate(Option_t *)
 501 {
 502   // Draw result to the screen
 503   // Called once at the end of the query
 504
 505 //  fHistPt = dynamic_cast<TH1F*> (GetOutputData(0));
 506 //  if (!fHistPt) {
 507 //    Printf("ERROR: fHistPt not available");
 508 //    return;
 509 //  }
 510
 511 //  TCanvas *c1 = new TCanvas("AliResonanceKinkPID","Pt MC",10,10,510,510);
 512 //  c1->cd(1)->SetLogy();
 513 //  fHistPt->DrawCopy("E");
 514 }
 515
 516 //____________________________________________________________________//
 517
 518 Float_t AliResonanceKinkPID::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 519   // Calculates the number of sigma to the vertex.
 520
 521   Float_t b[2];
 522   Float_t bRes[2];
 523   Float_t bCov[3];
 524
 525     esdTrack->GetImpactParameters(b,bCov);
 526
 527   if (bCov[0]<=0 || bCov[2]<=0) {
 528     //AliDebug(1, "Estimated b resolution lower or equal zero!");
 529     bCov[0]=0; bCov[2]=0;
 530   }
 531   bRes[0] = TMath::Sqrt(bCov[0]);
 532   bRes[1] = TMath::Sqrt(bCov[2]);
 533
 534   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 535
 536   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 537
 538   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 539
 540   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 541
 542   return d;
 543 }
 544
 545 //____________________________________________________________________//
 546
 547 const AliESDVertex* AliResonanceKinkPID::GetEventVertex(const AliESDEvent* esd) const
 548 {
 549   // Get the vertex
 550
 551   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 552
 553   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 554   else
 555   {
 556      vertex = esd->GetPrimaryVertexSPD();
 557      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 558      else
 559      return 0;
 560   }
 561 }
 562