]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG2/KINK/AliResonanceKink.cxx
From Stefano:
[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKink.cxx
index 1acfd9fab1995d2dc437cee8b60eebae57e7e8ca..5e95b2f7d7c3ce61f08896470af3a3a81c716537 100644 (file)
 //        Also, depending on the analysis mode (ESD or MC), fAnalysisType in the constructor must also be changed 
 //-----------------------------------------------------------------------------------------------------------------
 
-#include "TChain.h"
-#include "TTree.h"
 #include "TH2D.h"
 #include "TParticle.h"
 #include "TDatabasePDG.h"
-#include "TParticlePDG.h"
 #include "TF1.h"
 #include "TList.h"
-#include "TString.h"
-#include "AliMCEventHandler.h"
+
 #include "AliMCEvent.h"
 #include "AliResonanceKink.h"
 #include "AliESDkink.h"
@@ -43,20 +39,18 @@ ClassImp(AliResonanceKink)
 
 //________________________________________________________________________
 AliResonanceKink::AliResonanceKink() 
-  : 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), 
-  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), 
+  : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fetabins(0), floweta(0), fupeta(0), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(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), 
 fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
-, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
+, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0), fminKinkRadius(0), fmaxKinkRadius(0), fminQt(0), fmaxQt(0), fptbins(0), flowpt(0), fupperpt(0), fmaxAbsEtaCut(0)
 {
   // Constructor
 }
 
 //________________________________________________________________________
-AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t daughter1, Int_t daughter2, Int_t resonancePDG) 
-  : 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), 
-  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), 
+AliResonanceKink::AliResonanceKink(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t netabins, Float_t nloweta, Float_t nupeta, Int_t nptbins, Float_t nlowpt, Float_t nupperpt, Int_t daughter1, Int_t daughter2, Int_t resonancePDG) 
+  : TObject(), fDebug(0), fListOfHistos(0), fOpeningAngle(0), fInvariantMass(0), fInvMassTrue(0), fPhiBothKinks(0), fetabins(netabins), floweta(nloweta), fupeta(nupeta), fRecPt(0), fRecEta(0), fRecEtaPt(0), fSimPt(0), fSimEta(0), fSimEtaPt(0), fSimPtKink(0), fSimEtaKink(0),  fSimEtaPtKink(0), f1(0), f2(0), fAnalysisType(), fvtxz(0), fNbins(nbins), fLowX(nlowx), fHighX(nhighx), fdaughter1pdg(daughter1), fdaughter2pdg(daughter2), fresonancePDGcode(resonancePDG), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), 
 fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0), fMaxCov9(0), fMaxCov14(0) //, fTPCrefitFlag(kFALSE)
-, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0)
+, fInvmassPt(0), fInvmassPtTrue(0), fMCInvmassPt(0), fMCInvmassPtTrue(0), fminKinkRadius(0), fmaxKinkRadius(0), fminQt(0), fmaxQt(0), fptbins(nptbins), flowpt(nlowpt), fupperpt(nupperpt), fmaxAbsEtaCut(0)
 {
    // Constructor
   
@@ -66,17 +60,15 @@ fMaxDCAxy(0), fMaxDCAzaxis(0), fMinTPCclusters(0), fMaxChi2PerTPCcluster(0), fMa
    fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
    fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
 
-   fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
-   fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
-   fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); 
-   fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
-   fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); 
-   fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
-   fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
-   fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
-   fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);                
-   fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);  
-   fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
+   fRecPt=new TH1D("fRecPt"," ", nptbins, nlowpt, nupperpt);
+   fRecEta=new TH1D("fRecEta"," ",netabins, nloweta, nupeta);
+   fRecEtaPt=new TH2D("fRecEtaPt"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta); 
+   fSimPt=new TH1D("fSimPt"," ", nptbins, nlowpt, nupperpt);
+   fSimEta=new TH1D("fSimEta"," ", netabins, nloweta, nupeta); 
+   fSimEtaPt=new TH2D("fSimEtaPt"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta);
+   fSimPtKink=new TH1D("fSimPtKink"," ", nptbins, nlowpt, nupperpt);
+   fSimEtaKink=new TH1D("fSimEtaKink"," ", netabins, nloweta, nupeta);
+   fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta);                
    
    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);
    f1->SetParameter(0,0.493677);
@@ -113,8 +105,6 @@ AliResonanceKink:: ~AliResonanceKink()
  if(fSimPtKink) delete fSimPtKink;
  if(fSimEtaKink) delete fSimEtaKink;
  if(fSimEtaPtKink) delete fSimEtaPtKink;
- if(fhdr) delete fhdr;
- if(fhdz) delete fhdz;   
  if(fvtxz) delete fvtxz;  
  if(fInvmassPt) delete fInvmassPt; 
  if(fInvmassPtTrue) delete fInvmassPtTrue; 
@@ -140,8 +130,6 @@ TList* AliResonanceKink::GetHistogramList()
   fListOfHistos->Add(fSimPtKink);    
   fListOfHistos->Add(fSimEtaKink);   
   fListOfHistos->Add(fSimEtaPtKink);                                                           
-  fListOfHistos->Add(fhdr);
-  fListOfHistos->Add(fhdz);
   fListOfHistos->Add(fvtxz);
   fListOfHistos->Add(fInvmassPt);
   fListOfHistos->Add(fInvmassPtTrue);
@@ -152,7 +140,7 @@ TList* AliResonanceKink::GetHistogramList()
 }
 
 //________________________________________________________________________
-void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx)
+void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t nhighx, Int_t netabins, Float_t nloweta, Float_t nupeta, Int_t nptbins, Float_t nlowpt, Float_t nupperpt)
 {
   //  Initialisation of the output histograms
   fNbins=nbins; 
@@ -165,17 +153,15 @@ void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t
   fInvMassTrue=new TH1D("fInvMassTrue"," ",fNbins,fLowX,fHighX);
   fPhiBothKinks=new TH1D("fPhiBothKinks"," ",fNbins,fLowX,fHighX);  // Applicable for phi(1020)
 
-  fRecPt=new TH1D("fRecPt"," ", 50,0.0,5.0);
-  fRecEta=new TH1D("fRecEta"," ", 44,-1.1,1.1);
-  fRecEtaPt=new TH2D("fRecEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1); 
-  fSimPt=new TH1D("fSimPt"," ", 50,0.0,5.0);
-  fSimEta=new TH1D("fSimEta"," ", 44,-1.1,1.1); 
-  fSimEtaPt=new TH2D("fSimEtaPt"," ", 50,0.0,5.0, 44,-1.1,1.1);
-  fSimPtKink=new TH1D("fSimPtKink"," ", 50,0.0,5.0);
-  fSimEtaKink=new TH1D("fSimEtaKink"," ", 44,-1.1,1.1);
-  fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", 50,0.0,5.0, 44,-1.1,1.1);                
-  fhdr=new TH1D("fhdr"," ", 100,0.0,5.0);  
-  fhdz=new TH1D("fhdz"," ", 100,0.0,5.0);
+  fRecPt=new TH1D("fRecPt"," ", nptbins, nlowpt, nupperpt);
+  fRecEta=new TH1D("fRecEta"," ", netabins, nloweta, nupeta);
+  fRecEtaPt=new TH2D("fRecEtaPt"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta); 
+  fSimPt=new TH1D("fSimPt"," ", nptbins, nlowpt, nupperpt);
+  fSimEta=new TH1D("fSimEta"," ",netabins, nloweta, nupeta); 
+  fSimEtaPt=new TH2D("fSimEtaPt"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta);
+  fSimPtKink=new TH1D("fSimPtKink"," ", nptbins, nlowpt, nupperpt);
+  fSimEtaKink=new TH1D("fSimEtaKink"," ", netabins, nloweta, nupeta);
+  fSimEtaPtKink=new TH2D("fSimEtaPtKink"," ", nptbins, nlowpt, nupperpt, netabins, nloweta, nupeta);                
    
   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);
   f1->SetParameter(0,0.493677);
@@ -189,46 +175,47 @@ void AliResonanceKink::InitOutputHistograms(Int_t nbins, Float_t nlowx, Float_t
    
   fvtxz=new TH1D("fvtxz"," ", 100,-20.0,20.0);
   
-  fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
-  fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);  
-  fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
-  fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX,100,0.0,10.0);
+  fInvmassPt=new TH2D("fInvmassPt"," ",fNbins,fLowX,fHighX, nptbins, nlowpt, nupperpt);
+  fInvmassPtTrue=new TH2D("fInvmassPtTrue"," ",fNbins,fLowX,fHighX, nptbins, nlowpt, nupperpt);  
+  fMCInvmassPt=new TH2D("fMCInvmassPt"," ",fNbins,fLowX,fHighX,nptbins, nlowpt, nupperpt);
+  fMCInvmassPtTrue=new TH2D("fMCInvmassPtTrue"," ",fNbins,fLowX,fHighX, nptbins, nlowpt, nupperpt);
   
 }
-  
+
 //________________________________________________________________________
 void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent) 
 {
   // Main loop
   // Called for each event
   Int_t resonancePDGcode, antiresonancePDGcode;
+  Double_t daughter1pdgMass, daughter2pdgMass;
   
-  if (fdaughter1pdg==kdaughterKaon)  {
+  if (fdaughter1pdg==kKPlus)  {
     resonancePDGcode=fresonancePDGcode;
     antiresonancePDGcode=-fresonancePDGcode;
+    daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
+    daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass(); 
   }
-  if (fdaughter1pdg!=kdaughterKaon)  {
+  
+  if (fdaughter1pdg!=kKPlus)  {
     resonancePDGcode=-fresonancePDGcode;
     antiresonancePDGcode=fresonancePDGcode;
-  }  
+    daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
+    daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();   
+  }  // to ensure that daughter1pdgMass has always the kaon mass
+  
   if (fdaughter1pdg==fdaughter2pdg)  {
     resonancePDGcode=fresonancePDGcode;
     antiresonancePDGcode=fresonancePDGcode;
   }  
-
-   Double_t daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
-   Double_t daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
-   
   if (!esd) {
-    Printf("ERROR: fESD not available");
+    Printf("ERROR: esd not available");
     return;
-  }  
-
-    if (!mcEvent) {
-    Printf("ERROR: mcEvent not available");
-    return;
-  }  
-
+  }
+      if (!mcEvent) {
+      return;
+     }  
+  
   AliStack* stack=mcEvent->Stack();
 
   if(fAnalysisType == "MC") {
@@ -264,14 +251,14 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
        AliMCParticle   *mcDaughter1 = 0;
        AliMCParticle   *mcDaughter2 = 0;       
                
-        if(fdaughter1pdg==kdaughterKaon) {  
+        if(fdaughter1pdg==kKPlus) {  
           daughterParticle1=stack->Particle(firstD);
           daughterParticle2=stack->Particle(lastD);
          mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(firstD);  
          mcDaughter2= (AliMCParticle*) mcEvent->GetTrack(lastD);         
         }
         else 
-           if(fdaughter2pdg==kdaughterKaon) {
+           if(fdaughter2pdg==kKPlus) {
               daughterParticle1=stack->Particle(lastD);
               daughterParticle2=stack->Particle(firstD); 
              mcDaughter1= (AliMCParticle*) mcEvent->GetTrack(lastD);
@@ -280,16 +267,15 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
             
        if(TMath::Abs(daughterParticle1->GetPdgCode())!=321) continue;
        
-       TParticle       *daughters1Daughter = 0;
-       TParticle       *daughters2Daughter = 0;       
-       Int_t   mcProcessDaughters1Daughter = -999;
-       Int_t   mcProcessDaughters2Daughter = -999;       
+       TParticle      * daughters1Daughter=0;
+       TParticle      * daughters2Daughter=0;       
+       Int_t mcProcessDaughters1Daughter = -999;
+       Int_t mcProcessDaughters2Daughter = -999;       
        AliMCParticle *mcDaughters1Daughter = 0;
        
        if(mcDaughter1->Charge()==0) continue;
        if(mcDaughter2->Charge()==0) continue;       //accept resonance decays in two charged daughters
-
-       Int_t nDecayKaonDaughter = -99; 
+       Int_t nDecayKaonDaughter=-99; 
        for(Int_t ia=0; ia<daughterParticle1->GetNDaughters(); ia++) {
        if(((daughterParticle1->GetFirstDaughter()+ia)>0)&&((daughterParticle1->GetFirstDaughter()+ia)<stack->GetNtrack())) {
           daughters1Daughter=stack->Particle(daughterParticle1->GetFirstDaughter()+ia);
@@ -301,7 +287,7 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
         }
        }
        
-       Int_t nProcessDaughter = -99; 
+       Int_t nProcessDaughter=-99; 
        for(Int_t ib=0; ib<daughterParticle2->GetNDaughters(); ib++) {
         if(((daughterParticle2->GetFirstDaughter()+ib)>0)&&((daughterParticle2->GetFirstDaughter()+ib)<stack->GetNtrack())) {
           daughters2Daughter=stack->Particle(daughterParticle2->GetFirstDaughter()+ib);
@@ -313,7 +299,7 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
         }
        }  
        
-       Int_t numberOfCharged = 0;
+       Int_t numberOfCharged=0;
        if((mcProcessDaughters1Daughter==4)&&(nDecayKaonDaughter>=0)) {
          for(Int_t ic=nDecayKaonDaughter; ic<=daughterParticle1->GetLastDaughter(); ic++) {
           if ((ic>=0)&&(ic<stack->GetNtrack())) mcDaughters1Daughter= dynamic_cast<AliMCParticle*>(mcEvent->GetTrack(ic));
@@ -322,17 +308,16 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
          }
        }
        
-         if ((particle->Pt()>0.25)&&(TMath::Abs(particle->Eta())<0.9)) {
-
+                if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
+         if ((particle->Pt()>fMinPtTrackCut)&&(TMath::Abs(particle->Eta())<fmaxAbsEtaCut)) {
+         
          fSimEta->Fill(particle->Eta());
         fSimEtaPt->Fill(particle->Pt(), particle->Eta());
         fSimPt->Fill(particle->Pt());
-         fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());
-        if(numberOfCharged>=2) continue; // leave out kaon decays to more than one charged daughter
-          
+         fMCInvmassPtTrue->Fill(particle->GetMass(), particle->Pt());     
 
-         if((daughterParticle1->Pt()>0.25)&&(TMath::Abs(daughterParticle1->Eta())<0.9)&&(daughterParticle2->Pt()>0.25)&&(TMath::Abs(daughterParticle2->Eta())<0.9)) {
-          if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>120.)&&(daughters1Daughter->R()<220.)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>120.)&&(nProcessDaughter>0)))) { //below are the findable
+         if((daughterParticle1->Pt()>fMinPtTrackCut)&&(TMath::Abs(daughterParticle1->Eta())<fmaxAbsEtaCut)&&(daughterParticle2->Pt()>fMinPtTrackCut)&&(TMath::Abs(daughterParticle2->Eta())<fmaxAbsEtaCut)) {
+          if((mcProcessDaughters1Daughter==4)&&(daughters1Daughter->R()>fminKinkRadius)&&(daughters1Daughter->R()<fmaxKinkRadius)&&( (nProcessDaughter<0)||((daughters2Daughter->R()>fminKinkRadius)&&(nProcessDaughter>0)))) { //below are the findable
             fSimPtKink->Fill(particle->Pt());
             fSimEtaKink->Fill(particle->Eta());
             fSimEtaPtKink->Fill(particle->Pt(), particle->Eta());
@@ -345,8 +330,7 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
           }   //for the particular resonance
   } //MC loop
   
-  } // end fAnalysisType==MC
-  else 
+  } // end fAnalysisType==MC 
   
   if(fAnalysisType == "ESD") {
   const AliESDVertex* vertex = GetEventVertex(esd);
@@ -389,10 +373,17 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
     Int_t mumpdgpos = mumpos->GetPdgCode();
     
     Int_t indexKinkPos=trackpos->GetKinkIndex(0);
+    
+    if(indexKinkPos>0) continue;
+    
     Bool_t posKaonKinkFlag=0;
-    if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
     
-    if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
+    if(indexKinkPos<0) {
+      posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
+    
+      if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
+      if(posKaonKinkFlag==0) continue;
+    }
     
     if(indexKinkPos==0) {
 
@@ -429,10 +420,17 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
         Int_t mumpdgneg = mumneg->GetPdgCode();
        
        Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
+       
+       if(indexKinkNeg>0) continue;
+       
        Bool_t negKaonKinkFlag=0;
-       if(indexKinkNeg<0) negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
        
-       if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
+       if(indexKinkNeg<0) {
+         negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
+       
+         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
+         if(negKaonKinkFlag==0) continue;
+       }
        
        if(indexKinkNeg==0)  {
  
@@ -445,23 +443,29 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
        
        Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
 
-       if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
-        p4comb=anp4pos;
-        p4comb+=p4neg;
-        if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
+        if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
+         p4comb=anp4pos;
+         p4comb+=p4neg;
+         if((p4comb.Vect().Pt()>fMinPtTrackCut)&&(TMath::Abs(anp4pos.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(p4neg.Vect().Eta())<fmaxAbsEtaCut)&&(p4comb.Vect().Eta()<fmaxAbsEtaCut)) {        
+           if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
+         }
        }
                
        if(negKaonKinkFlag==1) {
          p4comb=p4pos;
           p4comb+=p4neg;
-         fInvariantMass->Fill(p4comb.M());
-         fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
-         if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
-          &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
-            fOpeningAngle->Fill(openingAngle);
-            fInvMassTrue->Fill(p4comb.M());
-           fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
-           if((TMath::Abs(p4pos.Vect().Eta())<1.1)&&(TMath::Abs(p4neg.Vect().Eta())<1.1)&&(p4comb.Vect().Eta()<1.1)) {
+         
+         if(p4comb.Vect().Pt()<=fMinPtTrackCut) continue;
+         
+         if((TMath::Abs(p4pos.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(p4neg.Vect().Eta())<fmaxAbsEtaCut)&&(p4comb.Vect().Eta()<fmaxAbsEtaCut)) {       
+           fInvariantMass->Fill(p4comb.M());
+           fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
+           if ((mumpdgpos==(antiresonancePDGcode))&&(mumpdgneg==(antiresonancePDGcode))&&(mumlabelpos==mumlabelneg)
+            &&(pdgpos==fdaughter2pdg)&&(pdgneg==(-fdaughter1pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)&&(mumlabelneg>=0)) {
+              fOpeningAngle->Fill(openingAngle);
+              fInvMassTrue->Fill(p4comb.M());
+             fInvmassPtTrue->Fill(p4comb.M(), p4comb.Vect().Pt());
+
              fRecPt->Fill(p4comb.Vect().Pt());
              fRecEta->Fill(p4comb.Vect().Eta());
              fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
@@ -475,17 +479,21 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
        if(posKaonKinkFlag==1) {
           anp4comb=anp4pos;
           anp4comb+=anp4neg;  
-         fInvariantMass->Fill(anp4comb.M());
-         fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
-         if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
-          &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
-            fOpeningAngle->Fill(openingAngle);
-            fInvMassTrue->Fill(p4comb.M());
-           fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
-            if((TMath::Abs(anp4neg.Vect().Eta())<1.1)&&(TMath::Abs(anp4pos.Vect().Eta())<1.1)&&(anp4comb.Vect().Eta()<1.1)) {  
-            fRecPt->Fill(anp4comb.Vect().Pt());
-            fRecEta->Fill(anp4comb.Vect().Eta());
-            fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
+         
+         if(anp4comb.Vect().Pt()<=fMinPtTrackCut) continue;      
+         
+          if((TMath::Abs(anp4neg.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(anp4pos.Vect().Eta())<fmaxAbsEtaCut)&&(anp4comb.Vect().Eta()<fmaxAbsEtaCut)) {        
+           fInvariantMass->Fill(anp4comb.M());
+           fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
+           if ((mumpdgpos==resonancePDGcode)&&(mumpdgneg==resonancePDGcode)&&(mumlabelpos==mumlabelneg)
+            &&(pdgpos==fdaughter1pdg)&&(pdgneg==(-fdaughter2pdg))&&(TMath::Abs(trackpos->GetLabel())>=0)&&(TMath::Abs(trackneg->GetLabel())>=0)&&(mumlabelpos>=0)  &&(mumlabelneg>=0)) {
+              fOpeningAngle->Fill(openingAngle);
+              fInvMassTrue->Fill(anp4comb.M());
+             fInvmassPtTrue->Fill(anp4comb.M(), anp4comb.Vect().Pt());
+       
+             fRecPt->Fill(anp4comb.Vect().Pt());
+             fRecEta->Fill(anp4comb.Vect().Eta());
+             fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta());
           }
 
          }
@@ -497,9 +505,170 @@ void AliResonanceKink::Analyse(AliESDEvent* esd, AliMCEvent* mcEvent)
   } //outer track loop 
   
   } // end fAnalysisType == ESD
-  
+    
 }      
 
+//____________________________________________________________________//
+void AliResonanceKink::Analyse(AliESDEvent* esd) 
+{
+  // Main loop
+  // Called for each event
+  Double_t daughter1pdgMass, daughter2pdgMass;
+  
+  if (fdaughter1pdg==kKPlus)  {
+    daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
+    daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass(); 
+  }
+  
+  if (fdaughter1pdg!=kKPlus)  {
+    daughter1pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
+    daughter2pdgMass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();   
+  }  // to ensure that daughter1pdgMass has always the kaon mass
+  
+  if (!esd) {
+    Printf("ERROR: esd not available");
+    return;
+  }
+ if(fAnalysisType == "DATA") {
+  const AliESDVertex* vertex = GetEventVertex(esd);
+  if(!vertex) return;
+  Double_t vtx[3];
+  vertex->GetXYZ(vtx);
+  fvtxz->Fill(vtx[2]);
+  Double_t ptrackpos[3], ptrackneg[3];
+  
+  TLorentzVector p4pos, anp4pos;
+  TLorentzVector p4neg, anp4neg;
+  TLorentzVector p4comb, anp4comb;
+  
+  for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
+    AliESDtrack* trackpos = esd->GetTrack(iTracks);
+    if (!trackpos) {
+      if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
+      continue;
+    }
+    if (trackpos->GetSign() < 0) continue;
+    
+    AliExternalTrackParam *tpcTrackpos = (AliExternalTrackParam *)trackpos->GetTPCInnerParam();
+    if(!tpcTrackpos) continue;
+    ptrackpos[0]=tpcTrackpos->Px();
+    ptrackpos[1]=tpcTrackpos->Py();   
+    ptrackpos[2]=tpcTrackpos->Pz();  
+    
+    Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
+    if(firstLevelAcceptPosTrack==kFALSE) continue;
+    
+    TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
+    
+    Int_t indexKinkPos=trackpos->GetKinkIndex(0);
+    
+    if(indexKinkPos>0) continue;
+    
+    Bool_t posKaonKinkFlag=0;
+    
+    if(indexKinkPos<0) {
+      posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
+    
+      if(posKaonKinkFlag==1) anp4pos.SetVectM(posTrackMom,daughter1pdgMass);
+      if(posKaonKinkFlag==0) continue;
+    }
+    
+    if(indexKinkPos==0) {
+
+    Bool_t secondLevelAcceptPosTrack=IsAcceptedForTrack(esd, vertex, trackpos);
+    if(secondLevelAcceptPosTrack==kFALSE) continue;
+
+      p4pos.SetVectM(posTrackMom, daughter2pdgMass);
+    
+    }
+       
+      for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
+        if(iTracks==j) continue;
+        AliESDtrack* trackneg=esd->GetTrack(j);
+        if (trackneg->GetSign() > 0) continue;
+       
+        AliExternalTrackParam *tpcTrackneg = (AliExternalTrackParam *)trackneg->GetTPCInnerParam();
+        if(!tpcTrackneg) continue;
+        ptrackneg[0]=tpcTrackneg->Px();
+        ptrackneg[1]=tpcTrackneg->Py();   
+        ptrackneg[2]=tpcTrackneg->Pz();  
+    
+        Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
+        if(firstLevelAcceptNegTrack==kFALSE) continue; 
+
+        TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
+       
+       Int_t indexKinkNeg=trackneg->GetKinkIndex(0);
+       
+       if(indexKinkNeg>0) continue;
+       
+       Bool_t negKaonKinkFlag=0;
+       
+       if(indexKinkNeg<0) {
+         negKaonKinkFlag=IsKink(esd, indexKinkNeg, negTrackMom);
+       
+         if(negKaonKinkFlag==1) p4neg.SetVectM(negTrackMom,daughter1pdgMass);
+         if(negKaonKinkFlag==0) continue;
+       }
+       
+       if(indexKinkNeg==0)  {
+          Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
+           if(secondLevelAcceptNegTrack==kFALSE) continue;  
+         
+         anp4neg.SetVectM(negTrackMom, daughter2pdgMass);
+       
+        }
+       
+       Double_t openingAngle=(ptrackpos[0]*ptrackneg[0]+ptrackpos[1]*ptrackneg[1]+ptrackpos[2]*ptrackneg[2])/(posTrackMom.Mag()*negTrackMom.Mag());
+
+        if((posKaonKinkFlag==1)&&(negKaonKinkFlag==1)) {
+         p4comb=anp4pos;
+         p4comb+=p4neg;
+         if((p4comb.Vect().Pt()>fMinPtTrackCut)&&(TMath::Abs(anp4pos.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(p4neg.Vect().Eta())<fmaxAbsEtaCut)&&(p4comb.Vect().Eta()<fmaxAbsEtaCut)) {        
+           if(openingAngle>0.6) fPhiBothKinks->Fill(p4comb.M());
+         }
+       }
+               
+       if(negKaonKinkFlag==1) {
+         p4comb=p4pos;
+          p4comb+=p4neg;
+         
+         if(p4comb.Vect().Pt()<=fMinPtTrackCut) continue;
+         
+         if((TMath::Abs(p4pos.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(p4neg.Vect().Eta())<fmaxAbsEtaCut)&&(p4comb.Vect().Eta()<fmaxAbsEtaCut)) {       
+           fInvariantMass->Fill(p4comb.M());
+           fInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
+           fRecPt->Fill(p4comb.Vect().Pt());
+           fRecEta->Fill(p4comb.Vect().Eta());
+           fRecEtaPt->Fill(p4comb.Vect().Perp(),p4comb.Vect().Eta());
+           }
+       }
+       
+       if(posKaonKinkFlag==1) {
+          anp4comb=anp4pos;
+          anp4comb+=anp4neg;  
+         
+         if(anp4comb.Vect().Pt()<=fMinPtTrackCut) continue;      
+         
+          if((TMath::Abs(anp4neg.Vect().Eta())<fmaxAbsEtaCut)&&(TMath::Abs(anp4pos.Vect().Eta())<fmaxAbsEtaCut)&&(anp4comb.Vect().Eta()<fmaxAbsEtaCut)) {        
+           fInvariantMass->Fill(anp4comb.M());
+           fInvmassPt->Fill(anp4comb.M(), anp4comb.Vect().Pt());
+           fRecPt->Fill(anp4comb.Vect().Pt());
+           fRecEta->Fill(anp4comb.Vect().Eta());
+           fRecEtaPt->Fill(anp4comb.Vect().Pt(), anp4comb.Vect().Eta()); 
+          }
+
+       }
+        
+      } //inner track loop 
+
+  } //outer track loop 
+  
+ } // end fAnalysisType == DATA
+
+}
 //____________________________________________________________________//
 Float_t AliResonanceKink::GetSigmaToVertex(AliESDtrack* esdTrack) const {
   // Calculates the number of sigma to the vertex.
@@ -585,7 +754,7 @@ const AliESDVertex* AliResonanceKink::GetEventVertex(const AliESDEvent* esd) con
       return kFALSE;
   }
   
-  if(gPt < fMinPtTrackCut) {
+  if(gPt <= fMinPtTrackCut) {
       if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
       return kFALSE;
   } 
@@ -630,7 +799,7 @@ Bool_t AliResonanceKink::IsAcceptedForTrack(AliESDEvent *localesd,
   } 
 
   if(nClustersTPC < fMinTPCclusters) {
-      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
+      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %d (min. requested: %d)", nClustersTPC, fMinTPCclusters);
       return kFALSE;
   } 
   
@@ -693,7 +862,7 @@ Bool_t AliResonanceKink::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3
          Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)*(p1XM-p2XM))+((p1YM-p2YM)*(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
          Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)*(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()*motherMfromKink.Mag());
 
-         if((kinkAngle>maxDecAngpimu)&&(qt>0.05)&&(qt<0.25)&&((kink->GetR()>110.)&&(kink->GetR()<230.))&&(TMath::Abs(trackMom.Eta())<1.1)&&(invariantMassKmu<0.6)) {
+         if((kinkAngle>maxDecAngpimu)&&(qt>fminQt)&&(qt<fmaxQt)&&((kink->GetR()>fminKinkRadius)&&(kink->GetR()<fmaxKinkRadius))&&(TMath::Abs(trackMom.Eta())<fmaxAbsEtaCut)&&(invariantMassKmu<0.6)) {
 
            if(trackMom.Mag()<=1.1) {
                return kTRUE;