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