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 "TDatabasePDG.h"
  28 #include "TParticlePDG.h"
  29 #include "TF1.h"
  30 #include "TList.h"
  31 #include "TString.h"
  32 #include "AliMCEventHandler.h"
  33 #include "AliMCEvent.h"
  34 #include "AliResonanceKink.h"
  35 #include "AliESDkink.h"
  36 #include "AliStack.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDEvent.h"
  39 #include "AliExternalTrackParam.h"
  40 #include "AliMCParticle.h"
  41
  42 ClassImp(AliResonanceKink)
  43
  44 //________________________________________________________________________
  45 AliResonanceKink::AliResonanceKink()
  46   : TObject(), fDebug(0), 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),
  47   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(0), fLowX(0), fHighX(0), fdaughter1pdg(0), fdaughter2pdg(0), fresonancePDGcode(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0),
  48 fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
  49 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
  50 {
  51   // Constructor
  52 }
  53
  54 //________________________________________________________________________
  55 AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG)
  56   : TObject(), fDebug(0), 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),
  57   fhdr(0), fhdz(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG), fMaxNSigmaToVertex(0), fMinPtTrackCut(0),
  58 fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
  59 , fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
  60 {
  61    // Constructor
  62
  63    fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
  64
  65    fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
  66    fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
  67    fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
  68
  69    fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
  70    fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
  71    fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
  72    fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
  73    fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
  74    fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
  75    fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
  76    fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
  77    fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
  78    fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
  79    fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
  80
  81    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);
  82    f1->SetParameter(0,0.493677);
  83    f1->SetParameter(1,0.9127037);
  84    f1->SetParameter(2,TMath::Pi());
  85
  86    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);
  87    f2->SetParameter(0,0.13957018);
  88    f2->SetParameter(1,0.2731374);
  89    f2->SetParameter(2,TMath::Pi());
  90
  91    fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
  92
  93    fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
  94    fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
  95    fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
  96    fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
  97 }
  98
  99 //________________________________________________________________________
 100 AliResonanceKink:: ~AliResonanceKink()
 101 {
 102  //  Destructor
 103  if(fOpeningAngle) delete fOpeningAngle;
 104  if(fInvariantMass) delete fInvariantMass;
 105  if(fInvMassTrue) delete fInvMassTrue;
 106  if(fPhiBothKinks) delete fPhiBothKinks;
 107  if(fRecPt) delete fRecPt;
 108  if(fRecEta) delete fRecEta;
 109  if(fRecEtaPt) delete fRecEtaPt;
 110  if(fSimPt) delete fSimPt;
 111  if(fSimEta) delete fSimEta;
 112  if(fSimEtaPt) delete fSimEtaPt;
 113  if(fSimPtKink) delete fSimPtKink;
 114  if(fSimEtaKink) delete fSimEtaKink;
 115  if(fSimEtaPtKink) delete fSimEtaPtKink;
 116  if(fhdr) delete fhdr;
 117  if(fhdz) delete fhdz;
 118  if(fvtxz) delete fvtxz;
 119  if(fInvmassPt) delete fInvmassPt;
 120  if(fInvmassPtTrue) delete fInvmassPtTrue;
 121  if(fMCInvmassPt) delete fMCInvmassPt;
 122  if(fMCInvmassPtTrue) delete fMCInvmassPtTrue;
 123 }
 124 //________________________________________________________________________
 125 TList* AliResonanceKink::GetHistogramList()
 126 {
 127   // Adding histograms to the list
 128   fListOfHistos=new TList();
 129
 130   fListOfHistos->Add(fOpeningAngle);
 131   fListOfHistos->Add(fInvariantMass);
 132   fListOfHistos->Add(fInvMassTrue);
 133   fListOfHistos->Add(fPhiBothKinks);
 134   fListOfHistos->Add(fRecPt);
 135   fListOfHistos->Add(fRecEta);
 136   fListOfHistos->Add(fRecEtaPt);
 137   fListOfHistos->Add(fSimPt);
 138   fListOfHistos->Add(fSimEta);
 139   fListOfHistos->Add(fSimEtaPt);
 140   fListOfHistos->Add(fSimPtKink);
 141   fListOfHistos->Add(fSimEtaKink);
 142   fListOfHistos->Add(fSimEtaPtKink);
 143   fListOfHistos->Add(fhdr);
 144   fListOfHistos->Add(fhdz);
 145   fListOfHistos->Add(fvtxz);
 146   fListOfHistos->Add(fInvmassPt);
 147   fListOfHistos->Add(fInvmassPtTrue);
 148   fListOfHistos->Add(fMCInvmassPt);
 149   fListOfHistos->Add(fMCInvmassPtTrue);
 150
 151   return fListOfHistos;
 152 }
 153
 154 //________________________________________________________________________
 155 void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
 156 {
 157   //  Initialisation of the output histograms
 158   fNbins=nbins;
 159   fLowX=nlowx;
 160   fHighX=nhighx;
 161
 162   fOpeningAngle=new TH1D("fOpeningAngle"," ", 100,-1.0,1.0);
 163
 164   fInvariantMass=new TH1D("fInvariantMass"," ",fNbins,fLowX,fHighX);
 165   fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
 166   fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
 167
 168   fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
 169   fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
 170   fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 171   fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
 172   fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1);
 173   fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
 174   fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
 175   fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
 176   fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);
 177   fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);
 178   fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
 179
 180   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);
 181   f1->SetParameter(0,0.493677);
 182   f1->SetParameter(1,0.9127037);
 183   f1->SetParameter(2,TMath::Pi());
 184
 185   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);
 186   f2->SetParameter(0,0.13957018);
 187   f2->SetParameter(1,0.2731374);
 188   f2->SetParameter(2,TMath::Pi());
 189
 190   fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
 191
 192   fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
 193   fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
 194   fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
 195   fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
 196
 197 }
 198
 199 //________________________________________________________________________
 200 void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
 201 {
 202   // Main loop
 203   // Called for each event
 204   Int_t resonancePDGcode, antiresonancePDGcode;
 205
 206   if (fdaughter1pdg==kdaughterKaon)  {
 207     resonancePDGcode=fresonancePDGcode;
 208     antiresonancePDGcode=-fresonancePDGcode;
 209   }
 210   if (fdaughter1pdg!=kdaughterKaon)  {
 211     resonancePDGcode=-fresonancePDGcode;
 212     antiresonancePDGcode=fresonancePDGcode;
 213   }
 214   if (fdaughter1pdg==fdaughter2pdg)  {
 215     resonancePDGcode=fresonancePDGcode;
 216     antiresonancePDGcode=fresonancePDGcode;
 217   }
 218
 219    Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
 220    Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
 221
 222   if (!esd) {
 223     Printf("ERROR: fESD not available");
 224     return;
 225   }
 226
 227     if (!mcEvent) {
 228     Printf("ERROR: mcEvent not available");
 229     return;
 230   }
 231
 232   AliStack* stack=mcEvent->Stack();
 233
 234   if(fAnalysisType == "MC") {
 235
 236   const AliESDVertex* vertex = GetEventVertex(esd);
 237   if(!vertex) return;
 238   Double_t vtx[3];
 239   vertex->GetXYZ(vtx);
 240   fvtxz->Fill(vtx[2]);
 241
 242   for (Int_t iMc = 0; iMc < stack->GetNprimary(); ++iMc)
 243   {
 244     TParticle* particle = stack->Particle(iMc);
 245
 246     if (!particle)
 247     {
 248       if (fDebug > 0) Printf("UNEXPECTED: particle with label %d not found in stack (mc loop)", iMc);
 249       continue;
 250     }
 251
 252        if(TMath::Abs(particle->GetPdgCode())==fresonancePDGcode) {
 253
 254        if(particle->GetNDaughters()>2) continue;
 255
 256        Int_t firstD=particle->GetFirstDaughter();
 257        Int_t lastD=particle->GetLastDaughter();
 258
 259        if((firstD<0)||(firstD>=stack->GetNtrack())) continue;
 260        if((lastD<0)||(lastD>=stack->GetNtrack())) continue;
 261
 262        TParticle *daughterParticle1;
 263        TParticle *daughterParticle2;
 264        AliMCParticle *mcDaughter1;
 265        AliMCParticle *mcDaughter2;
 266
 267         if(fdaughter1pdg==kdaughterKaon) {
 268           daughterParticle1=stack->Particle(firstD);
 269           daughterParticle2=stack->Particle(lastD);
 270           mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);
 271           mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);
 272         }
 273         else {
 274             if(fdaughter2pdg==kdaughterKaon) {
 275               daughterParticle1=stack->Particle(lastD);
 276               daughterParticle2=stack->Particle(firstD);
 277               mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD);
 278               mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(firstD);
 279              }    //to ensure that the first daughter is always the kaon
 280             else continue;
 281         }
 282
 283        TParticle* daughters1Daughter=0;
 284        TParticle* daughters2Daughter=0;
 285        Int_t mcProcessDaughters1Daughter = -999;
 286        Int_t mcProcessDaughters2Daughter = -999;
 287        AliMCParticle *mcDaughters1Daughter;
 288
 289        if(mcDaughter1->Charge()==0) continue;
 290        if(mcDaughter2->Charge()==0) continue;       //accept resonance decays in two charged daughters
 291
 292        Int_t nDecayKaonDaughter=-99;
 293        for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) {
 294        if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) {
 295            daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
 296            mcProcessDaughters1Daughter=daughters1Daughter->GetUniqueID();
 297            if(mcProcessDaughters1Daughter==4) {
 298              nDecayKaonDaughter=daughterParticle1->GetFirstDaughter()+ia;
 299              break;
 300            }
 301         }
 302        }
 303
 304        Int_t nProcessDaughter=-99;
 305        for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) {
 306         if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) {
 307            daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib);
 308            mcProcessDaughters2Daughter=daughters2Daughter->GetUniqueID();
 309            if((mcProcessDaughters2Daughter==4)||(mcProcessDaughters2Daughter==13)) {
 310              nProcessDaughter=mcProcessDaughters2Daughter;
 311              break;
 312            }
 313         }
 314        }
 315
 316        Int_t numberOfCharged=0;
 317        if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) {
 318          for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) {
 319            if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic));
 320             else continue;
 321             if(mcDaughters1Daughter->Charge()!=0) numberOfCharged=numberOfCharged+1;
 322          }
 323        }
 324
 325          if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
 326
 327          fSimEta->Fill(particle->Eta());
 328          fSimEtaPt->Fill(particle->Pt(), particle->Eta());
 329          fSimPt->Fill(particle->Pt());
 330          fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());
 331          if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
 332
 333
 334          if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) {
 335            if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable
 336              fSimPtKink->Fill(particle->Pt());
 337              fSimEtaKink->Fill(particle->Eta());
 338              fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
 339              fMCInvmassPt->Fill(particle->GetMass(), particle->Pt());
 340            }
 341          }
 342
 343        } // for the generated spectra
 344
 345           }   //for the particular resonance
 346   } //MC loop
 347
 348   } // end fAnalysisType==MC
 349   else
 350
 351   if(fAnalysisType == "ESD") {
 352   const AliESDVertex* vertex = GetEventVertex(esd);
 353   if(!vertex) return;
 354   Double_t vtx[3];
 355   vertex->GetXYZ(vtx);
 356   fvtxz->Fill(vtx[2]);
 357   Double_t ptrackpos[3], ptrackneg[3];
 358
 359   TLorentzVector p4pos, anp4pos;
 360   TLorentzVector p4neg, anp4neg;
 361   TLorentzVector p4comb, anp4comb;
 362
 363   for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
 364     AliESDtrack* trackpos = esd->GetTrack(iTracks);
 365     if (!trackpos) {
 366       if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
 367       continue;
 368     }
 369     if (trackpos->GetSign() < 0) continue;
 370
 371     AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
 372     if(!tpcTrackpos) continue;
 373     ptrackpos[0]=tpcTrackpos->Px();
 374     ptrackpos[1]=tpcTrackpos->Py();
 375     ptrackpos[2]=tpcTrackpos->Pz();
 376
 377     Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
 378     if(firstLevelAcceptPosTrack==kFALSE) continue;
 379
 380     TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
 381
 382     TParticle * partpos = stack->Particle(TMath::Abs(trackpos->GetLabel()));
 383     if (!partpos) continue;
 384     Int_t pdgpos = partpos->GetPdgCode();
 385     Int_t mumlabelpos=partpos->GetFirstMother();
 386     mumlabelpos = TMath::Abs(mumlabelpos);
 387     TParticle * mumpos=stack->Particle(mumlabelpos);
 388     if (!mumpos) continue;
 389     Int_t mumpdgpos = mumpos->GetPdgCode();
 390
 391     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
 392     Bool_t posKaonKinkFlag=0;
 393     if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
 394
 395     if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
 396
 397     if(indexKinkPos==0) {
 398
 399     Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
 400     if(secondLevelAcceptPosTrack==kFALSE) continue;
 401
 402       p4pos.SetVectM(posTrackMom, daughter2pdgMass);
 403
 404     }
 405
 406       for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
 407         if(iTracks==j) continue;
 408         AliESDtrack* trackneg=esd->GetTrack(j);
 409         if (trackneg->GetSign() > 0) continue;
 410
 411         AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
 412         if(!tpcTrackneg) continue;
 413         ptrackneg[0]=tpcTrackneg->Px();
 414         ptrackneg[1]=tpcTrackneg->Py();
 415         ptrackneg[2]=tpcTrackneg->Pz();
 416
 417         Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
 418         if(firstLevelAcceptNegTrack==kFALSE) continue;
 419
 420         TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
 421
 422         TParticle * partneg = stack->Particle(TMath::Abs(trackneg->GetLabel()));
 423         if (!partneg) continue;
 424         Int_t pdgneg = partneg->GetPdgCode();
 425         Int_t mumlabelneg=partneg->GetFirstMother();
 426         mumlabelneg = TMath::Abs(mumlabelneg);
 427         TParticle * mumneg=stack->Particle(mumlabelneg);
 428         if (!mumneg) continue;
 429         Int_t mumpdgneg = mumneg->GetPdgCode();
 430
 431         Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
 432         Bool_t negKaonKinkFlag=0;
 433         if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
 434
 435         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
 436
 437         if(indexKinkNeg==0)  {
 438
 439            Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
 440            if(secondLevelAcceptNegTrack==kFALSE) continue;
 441
 442           anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
 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((posKaonKinkFlag==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           fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
 459           if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
 460           &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
 461             fOpeningAngle->Fill(openingAngle);
 462             fInvMassTrue->Fill(p4comb.M());
 463             fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
 464             if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
 465               fRecPt->Fill(p4comb.Vect().Pt());
 466               fRecEta->Fill(p4comb.Vect().Eta());
 467               fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
 468
 469             }
 470
 471            }
 472
 473         }
 474
 475         if(posKaonKinkFlag==1) {
 476           anp4comb=anp4pos;
 477           anp4comb+=anp4neg;
 478           fInvariantMass->Fill(anp4comb.M());
 479           fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
 480           if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
 481           &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
 482             fOpeningAngle->Fill(openingAngle);
 483             fInvMassTrue->Fill(p4comb.M());
 484             fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
 485             if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {
 486              fRecPt->Fill(anp4comb.Vect().Pt());
 487              fRecEta->Fill(anp4comb.Vect().Eta());
 488              fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
 489            }
 490
 491          }
 492
 493         }
 494
 495       } //inner track loop
 496
 497   } //outer track loop
 498
 499   } // end fAnalysisType == ESD
 500
 501 }
 502
 503 //____________________________________________________________________//
 504 Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
 505   // Calculates the number of sigma to the vertex.
 506
 507   Float_t b[2];
 508   Float_t bRes[2];
 509   Float_t bCov[3];
 510
 511     esdTrack->GetImpactParametersTPC(b,bCov);
 512
 513   if (bCov[0]<=0 || bCov[2]<=0) {
 514     bCov[0]=0; bCov[2]=0;
 515   }
 516
 517   bRes[0] = TMath::Sqrt(bCov[0]);
 518   bRes[1] = TMath::Sqrt(bCov[2]);
 519
 520   if (bRes[0] == 0 || bRes[1] ==0) return -1;
 521
 522   Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
 523
 524   if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
 525
 526   d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
 527
 528   return d;
 529 }
 530
 531 //________________________________________________________________________
 532 const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) const
 533 {
 534   // Get the vertex
 535
 536   const AliESDVertex* vertex = esd->GetPrimaryVertex();
 537
 538   if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 539   else
 540   {
 541      vertex = esd->GetPrimaryVertexSPD();
 542      if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
 543      else
 544      return 0;
 545   }
 546 }
 547
 548 //________________________________________________________________________
 549
 550  Bool_t AliResonanceKink::IsAcceptedForKink(AliESDEvent *localesd,
 551             const AliESDVertex *localvertex, AliESDtrack* localtrack) {
 552    // Apply the selections for each kink
 553
 554   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 555   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 556   Double_t dca3D = 0.0;
 557
 558   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
 559   if(!tpcTrack) {
 560     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 561     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 562     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 563   }
 564   else {
 565     gPt = tpcTrack->Pt();
 566     gPx = tpcTrack->Px();
 567     gPy = tpcTrack->Py();
 568     gPz = tpcTrack->Pz();
 569     tpcTrack->PropagateToDCA(localvertex,
 570                localesd->GetMagneticField(),100.,dca,cov);
 571   }
 572
 573   if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
 574       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",   GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
 575       return kFALSE;
 576   }
 577
 578   if(TMath::Abs(dca[0]) > fMaxDCAxy) {
 579       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(xy) (TPC) of %lf (max. requested: %lf)", TMath::Abs(dca[0]), fMaxDCAxy);
 580       return kFALSE;
 581   }
 582
 583   if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
 584       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a value of dca(z) of %lf (max. requested: %lf)", TMath::Abs(dca[1]), fMaxDCAzaxis);
 585       return kFALSE;
 586   }
 587
 588   if(gPt < fMinPtTrackCut) {
 589       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
 590       return kFALSE;
 591   }
 592
 593   return kTRUE;
 594 }
 595
 596 //________________________________________________________________________
 597 Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd,                                                                                                                                          const AliESDVertex *localvertex, AliESDtrack *localtrack) {
 598    // Apply the selections for each track
 599
 600   Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
 601   Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
 602   Double_t dca3D = 0.0;
 603
 604   AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
 605   if(!tpcTrack) {
 606     gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
 607     dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
 608     cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
 609   }
 610   else {
 611     gPt = tpcTrack->Pt();
 612     gPx = tpcTrack->Px();
 613     gPy = tpcTrack->Py();
 614     gPz = tpcTrack->Pz();
 615     tpcTrack->PropagateToDCA(localvertex,
 616                localesd->GetMagneticField(),100.,dca,cov);
 617   }
 618
 619   Int_t fcls[200];
 620   Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
 621   Float_t chi2perTPCcluster=-1.0;
 622   if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
 623
 624   Double_t extCov[15];
 625   localtrack->GetExternalCovariance(extCov);
 626
 627   if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
 628       if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
 629       return kFALSE;
 630   }
 631
 632   if(nClustersTPC < fMinTPCclusters) {
 633       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
 634       return kFALSE;
 635   }
 636
 637   if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
 638       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
 639       return kFALSE;
 640   }
 641
 642   if(extCov[0] > fMaxCov0) {
 643       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
 644       return kFALSE;
 645   }
 646
 647   if(extCov[2] > fMaxCov2) {
 648       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
 649       return kFALSE;
 650   }
 651
 652   if(extCov[5] > fMaxCov5) {
 653       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
 654       return kFALSE;
 655   }
 656
 657   if(extCov[9] > fMaxCov9) {
 658       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
 659       return kFALSE;
 660   }
 661
 662   if(extCov[14] > fMaxCov14) {
 663       if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
 664       return kFALSE;
 665   }
 666
 667   return kTRUE;
 668
 669 }
 670
 671 //________________________________________________________________________
 672 Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
 673 {
 674    // Test some kinematical criteria for each kink
 675
 676          AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
 677          const TVector3 motherMfromKink(kink->GetMotherP());
 678          const TVector3 daughterMKink(kink->GetDaughterP());
 679          Float_t qt=kink->GetQt();
 680
 681          Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
 682          Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
 683
 684          Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
 685
 686          Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
 687          Float_t p1XM= motherMfromKink.Px();
 688          Float_t p1YM= motherMfromKink.Py();
 689          Float_t p1ZM= motherMfromKink.Pz();
 690          Float_t p2XM= daughterMKink.Px();
 691          Float_t p2YM= daughterMKink.Py();
 692          Float_t p2ZM= daughterMKink.Pz();
 693          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
 694          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 695
 696          if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
 697
 698            if(trackMom.Mag()<=1.1) {
 699                 return kTRUE;
 700            }
 701            else
 702            if (kinkAngle<maxDecAngKmu) {
 703                 return kTRUE;
 704            }
 705          }
 706          return kFALSE;
 707 }