[u/mrichter/AliRoot.git] / PWG2 / KINK / AliResonanceKinkLikeSign.cxx

/**************************************************************************
 * Author: Paraskevi Ganoti, University of Athens (pganoti@phys.uoa.gr)   *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//----------------------------------------------------------------------------------------------------------------------
//                        class AliResonanceKinkLikeSign
//        Example of an analysis task for producing a like-sign background for resonances having at least one 
//        kaon-kink in their decay products. 
//        Background is calculated from a positive kaon kink and a positive track but other possibilities are feasible.
//----------------------------------------------------------------------------------------------------------------------

#include "AliESDEvent.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TF1.h"
#include "TDatabasePDG.h"
#include "TLorentzVector.h"
#include "AliESDkink.h"

#include "AliResonanceKinkLikeSign.h"

ClassImp(AliResonanceKinkLikeSign)

//________________________________________________________________________
AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name) 
  : AliAnalysisTaskSE(name), fDebug(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0), 
fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0), fdaughter1pdg(0), fdaughter2pdg(0), fnbins(0), fnlowx(0), fnhighx(0), floweta(0), fuppereta(0), fminKinkRadius(0), fmaxKinkRadius(0), fminQt(0), fmaxQt(0), fptbins(0), flowpt(0), fupperpt(0)

{
  // Constructor

  // Define input and output slots here
  // Input slot #0 works with a TChain
  DefineOutput(1, TList::Class());
}

//________________________________________________________________________
void AliResonanceKinkLikeSign::UserCreateOutputObjects() 
{
  // Create histograms
  // Called once
  
   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);
   f1->SetParameter(1,0.9127037);
   f1->SetParameter(2,TMath::Pi());

   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);
   f2->SetParameter(0,0.13957018);
   f2->SetParameter(1,0.2731374);
   f2->SetParameter(2,TMath::Pi());
  
   fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", fnbins, fnlowx, fnhighx);
   fLikeSignInvmassPt=new TH2D("fLikeSignInvmassPt"," ", fnbins, fnlowx, fnhighx, fptbins, flowpt, fupperpt);
 
   fListOfHistos=new TList();
   fListOfHistos->Add(fPosKaonLikeSign);
   fListOfHistos->Add(fLikeSignInvmassPt);     

}

//________________________________________________________________________
void AliResonanceKinkLikeSign::UserExec(Option_t *) 
{
  // Main loop
  // Called for each event
 
  AliVEvent *event = InputEvent();
  if (!event) {
     Printf("ERROR: Could not retrieve event");
     return;
  }

  AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
  if (!esd) {
     Printf("ERROR: Could not retrieve esd");
     return;
  }

   Double_t daughter1Mass, daughter2Mass;
  
   if (fdaughter1pdg==kKPlus)  {
     daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
     daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass(); 
   }
  
   if (fdaughter1pdg!=kKPlus)  {
     daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
     daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();   
   }  // to ensure that daughter1Mass has always the kaon mass

  if (!esd) {
    Printf("ERROR: esd not available");
    return;
  }  

  const AliESDVertex* vertex = GetEventVertex(esd);
  if(!vertex) return;
  
  Double_t ptrackpos[3], ptrackneg[3];
  
  TLorentzVector p4pos;
  TLorentzVector p4neg;
  TLorentzVector p4comb; 

  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);
    Bool_t posKaonKinkFlag=0;
    if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
    if(posKaonKinkFlag==0) continue; 
    if(posKaonKinkFlag==1) p4pos.SetVectM(posTrackMom, daughter1Mass);

    for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
        if(iTracks==j) continue;
        AliESDtrack* trackneg=esd->GetTrack(j);
        if (!trackneg) {
          if (fDebug > 0) Printf("ERROR: Could not receive track %d", j);
          continue;
        } 
        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]);

	Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
        if(secondLevelAcceptNegTrack==kFALSE) continue;  
	 	
      	p4neg.SetVectM(negTrackMom, daughter2Mass);  
	
	p4comb=p4pos;
	p4comb+=p4neg;
	
	if(p4comb.Vect().Pt()<=fMinPtTrackCut) continue;	
	
        if((TMath::Abs(p4pos.Vect().Eta())<fuppereta)&&(TMath::Abs(p4neg.Vect().Eta())<fuppereta)&&(p4comb.Vect().Eta()<fuppereta)) {
	
	  fPosKaonLikeSign->Fill(p4comb.M());
	  fLikeSignInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
	
	}
	
      } //inner track loop 

  } //outer track loop 
    
    PostData(1, fListOfHistos);
}      

//________________________________________________________________________
void AliResonanceKinkLikeSign::Terminate(Option_t *) 
{
 
}

//____________________________________________________________________//

Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
  // Calculates the number of sigma to the vertex.
  
  Float_t b[2];
  Float_t bRes[2];
  Float_t bCov[3];

    esdTrack->GetImpactParameters(b,bCov);
  
  if (bCov[0]<=0 || bCov[2]<=0) {
    //AliDebug(1, "Estimated b resolution lower or equal zero!");
    bCov[0]=0; bCov[2]=0;
  }
  bRes[0] = TMath::Sqrt(bCov[0]);
  bRes[1] = TMath::Sqrt(bCov[2]);
  
  if (bRes[0] == 0 || bRes[1] ==0) return -1;
  
  Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
  
  if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
  
  d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
  
  return d;
}

//____________________________________________________________________//

const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const

{
  // Get the vertex

  const AliESDVertex* vertex = esd->GetPrimaryVertex();

  if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
  else
  {
     vertex = esd->GetPrimaryVertexSPD();
     if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
     else
     return 0;
  }
}

//________________________________________________________________________

 Bool_t AliResonanceKinkLikeSign::IsAcceptedForKink(AliESDEvent *localesd,
            const AliESDVertex *localvertex, AliESDtrack* localtrack) {
   // Apply the selections for each kink

  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
  Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
  Double_t dca3D = 0.0;
  
  AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
  if(!tpcTrack) {
    gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
    dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
    cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
  }
  else {
    gPt = tpcTrack->Pt();
    gPx = tpcTrack->Px();
    gPy = tpcTrack->Py();
    gPz = tpcTrack->Pz();
    tpcTrack->PropagateToDCA(localvertex,
    	       localesd->GetMagneticField(),100.,dca,cov);
  }
  
  if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
      if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)",   GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
      return kFALSE;
  }
  
  if(TMath::Abs(dca[0]) > fMaxDCAxy) {
      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);
      return kFALSE;
  }
    
  if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
      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);
      return kFALSE;
  }
  
  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;
  } 
  
  return kTRUE;
}

//________________________________________________________________________
Bool_t AliResonanceKinkLikeSign::IsAcceptedForTrack(AliESDEvent *localesd,                                                                                                                                          const AliESDVertex *localvertex, AliESDtrack *localtrack) {
   // Apply the selections for each track

  Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
  Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};  //The impact parameters and their covariance.
  Double_t dca3D = 0.0;
  
  AliExternalTrackParam *tpcTrack = (AliExternalTrackParam *)localtrack->GetTPCInnerParam();
  if(!tpcTrack) {
    gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
    dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
    cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
  }
  else {
    gPt = tpcTrack->Pt();
    gPx = tpcTrack->Px();
    gPy = tpcTrack->Py();
    gPz = tpcTrack->Pz();
    tpcTrack->PropagateToDCA(localvertex,
    	       localesd->GetMagneticField(),100.,dca,cov);
  }
  
  Int_t fcls[200];
  Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
  Float_t chi2perTPCcluster=-1.0;
  if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
  
  Double_t extCov[15];
  localtrack->GetExternalCovariance(extCov);
  
  if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
      return kFALSE;
  } 

  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);
      return kFALSE;
  } 
  
  if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
      return kFALSE;
  } 

  if(extCov[0] > fMaxCov0) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
      return kFALSE;
  }
  
  if(extCov[2] > fMaxCov2) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
      return kFALSE;
  }
    
  if(extCov[5] > fMaxCov5) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
      return kFALSE;
  }  
  
  if(extCov[9] > fMaxCov9) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
      return kFALSE;
  }  
  
  if(extCov[14] > fMaxCov14) {
      if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
      return kFALSE;
  } 
 
  return kTRUE;

}

//________________________________________________________________________
Bool_t AliResonanceKinkLikeSign::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom) 
{
   // Test some kinematical criteria for each kink

	 AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
	 const TVector3 motherMfromKink(kink->GetMotherP());
	 const TVector3 daughterMKink(kink->GetDaughterP());
	 Float_t qt=kink->GetQt();

	 Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
	 Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);

         Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
	 
	 Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()*daughterMKink.Mag()+0.105658*0.105658);
	 Float_t p1XM= motherMfromKink.Px();
         Float_t p1YM= motherMfromKink.Py();
         Float_t p1ZM= motherMfromKink.Pz();
         Float_t p2XM= daughterMKink.Px();
         Float_t p2YM= daughterMKink.Py();
         Float_t p2ZM= daughterMKink.Pz();
         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>fminQt)&&(qt<fmaxQt)&&((kink->GetR()>fminKinkRadius)&&(kink->GetR()<fmaxKinkRadius))&&(TMath::Abs(trackMom.Eta())<fuppereta)&&(invariantMassKmu<0.6)) {

           if(trackMom.Mag()<=1.1) {
		return kTRUE;
           }
	   else 
	   if (kinkAngle<maxDecAngKmu) {
		return kTRUE;
	   }
	 }
	 return kFALSE;
}
Commit	Line	Data
10eaad41	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
894840ad	14	//----------------------------------------------------------------------------------------------------------------------
10eaad41	15	// class AliResonanceKinkLikeSign
	16	// Example of an analysis task for producing a like-sign background for resonances having at least one
	17	// kaon-kink in their decay products.
894840ad	18	// Background is calculated from a positive kaon kink and a positive track but other possibilities are feasible.
894840ad	19	//----------------------------------------------------------------------------------------------------------------------
10eaad41	20
be1a7181	21	#include "AliESDEvent.h"
10eaad41	22	#include "TH1D.h"
7ccf0419	23	#include "TH2D.h"
be1a7181	24	#include "TF1.h"
be1a7181	25	#include "TDatabasePDG.h"
10eaad41	26	#include "TLorentzVector.h"
be1a7181	27	#include "AliESDkink.h"
10eaad41	28
10eaad41	29	#include "AliResonanceKinkLikeSign.h"
10eaad41	30
	31	ClassImp(AliResonanceKinkLikeSign)
	32
10eaad41	33	//________________________________________________________________________
10eaad41	34	AliResonanceKinkLikeSign::AliResonanceKinkLikeSign(const char *name)
be1a7181	35	: AliAnalysisTaskSE(name), fDebug(0), fListOfHistos(0), f1(0), f2(0), fPosKaonLikeSign(0), fLikeSignInvmassPt(0), fMaxNSigmaToVertex(0), fMinPtTrackCut(0), fMaxDCAxy(0), fMaxDCAzaxis(0),
894840ad	36	fMinTPCclusters(0),fMaxChi2PerTPCcluster(0), fMaxCov0(0), fMaxCov2(0), fMaxCov5(0) , fMaxCov9(0), fMaxCov14(0), fdaughter1pdg(0), fdaughter2pdg(0), fnbins(0), fnlowx(0), fnhighx(0), floweta(0), fuppereta(0), fminKinkRadius(0), fmaxKinkRadius(0), fminQt(0), fmaxQt(0), fptbins(0), flowpt(0), fupperpt(0)
10eaad41	37
	38	{
	39	// Constructor
	40
	41	// Define input and output slots here
	42	// Input slot #0 works with a TChain
10eaad41	43	DefineOutput(1, TList::Class());
	44	}
	45
	46	//________________________________________________________________________
be1a7181	47	void AliResonanceKinkLikeSign::UserCreateOutputObjects()
10eaad41	48	{
	49	// Create histograms
	50	// Called once
	51
	52	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);
	53	f1->SetParameter(0,0.493677);
	54	f1->SetParameter(1,0.9127037);
	55	f1->SetParameter(2,TMath::Pi());
	56
	57	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);
	58	f2->SetParameter(0,0.13957018);
	59	f2->SetParameter(1,0.2731374);
	60	f2->SetParameter(2,TMath::Pi());
	61
be1a7181	62	fPosKaonLikeSign=new TH1D("fPosKaonLikeSign"," ", fnbins, fnlowx, fnhighx);
894840ad	63	fLikeSignInvmassPt=new TH2D("fLikeSignInvmassPt"," ", fnbins, fnlowx, fnhighx, fptbins, flowpt, fupperpt);
be1a7181	64
10eaad41	65	fListOfHistos=new TList();
7ccf0419	66	fListOfHistos->Add(fPosKaonLikeSign);
7ccf0419	67	fListOfHistos->Add(fLikeSignInvmassPt);
10eaad41	68
	69	}
	70
	71	//________________________________________________________________________
be1a7181	72	void AliResonanceKinkLikeSign::UserExec(Option_t *)
10eaad41	73	{
	74	// Main loop
	75	// Called for each event
	76
be1a7181	77	AliVEvent *event = InputEvent();
	78	if (!event) {
	79	Printf("ERROR: Could not retrieve event");
	80	return;
	81	}
10eaad41	82
be1a7181	83	AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event);
	84	if (!esd) {
	85	Printf("ERROR: Could not retrieve esd");
	86	return;
	87	}
f2a16059	88
be1a7181	89	Double_t daughter1Mass, daughter2Mass;
	90
	91	if (fdaughter1pdg==kKPlus) {
	92	daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
	93	daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
	94	}
	95
	96	if (fdaughter1pdg!=kKPlus) {
	97	daughter1Mass=TDatabasePDG::Instance()->GetParticle(fdaughter2pdg)->Mass();
	98	daughter2Mass=TDatabasePDG::Instance()->GetParticle(fdaughter1pdg)->Mass();
	99	} // to ensure that daughter1Mass has always the kaon mass
	100
	101	if (!esd) {
	102	Printf("ERROR: esd not available");
10eaad41	103	return;
	104	}
	105
be1a7181	106	const AliESDVertex* vertex = GetEventVertex(esd);
10eaad41	107	if(!vertex) return;
	108
	109	Double_t ptrackpos[3], ptrackneg[3];
	110
	111	TLorentzVector p4pos;
	112	TLorentzVector p4neg;
	113	TLorentzVector p4comb;
	114
be1a7181	115	for (Int_t iTracks = 0; iTracks < esd->GetNumberOfTracks(); iTracks++) {
	116
	117	AliESDtrack* trackpos = esd->GetTrack(iTracks);
10eaad41	118	if (!trackpos) {
1ac5c71d	119	if (fDebug > 0) Printf("ERROR: Could not receive track %d", iTracks);
10eaad41	120	continue;
10eaad41	121	}
7ccf0419	122	if (trackpos->GetSign() < 0) continue;
10eaad41	123
be1a7181	124	AliExternalTrackParam tpcTrackpos = (AliExternalTrackParam )trackpos->GetTPCInnerParam();
	125	if(!tpcTrackpos) continue;
	126	ptrackpos[0]=tpcTrackpos->Px();
	127	ptrackpos[1]=tpcTrackpos->Py();
	128	ptrackpos[2]=tpcTrackpos->Pz();
10eaad41	129
be1a7181	130	Bool_t firstLevelAcceptPosTrack=IsAcceptedForKink(esd, vertex, trackpos);
be1a7181	131	if(firstLevelAcceptPosTrack==kFALSE) continue;
10eaad41	132
10eaad41	133	TVector3 posTrackMom(ptrackpos[0],ptrackpos[1],ptrackpos[2]);
10eaad41	134
be1a7181	135	Int_t indexKinkPos=trackpos->GetKinkIndex(0);
	136	Bool_t posKaonKinkFlag=0;
	137	if(indexKinkPos<0) posKaonKinkFlag=IsKink(esd, indexKinkPos, posTrackMom);
	138	if(posKaonKinkFlag==0) continue;
	139	if(posKaonKinkFlag==1) p4pos.SetVectM(posTrackMom, daughter1Mass);
7ccf0419	140
be1a7181	141	for (Int_t j=0; j<esd->GetNumberOfTracks(); j++) {
10eaad41	142	if(iTracks==j) continue;
be1a7181	143	AliESDtrack* trackneg=esd->GetTrack(j);
	144	if (!trackneg) {
	145	if (fDebug > 0) Printf("ERROR: Could not receive track %d", j);
	146	continue;
	147	}
7ccf0419	148	if (trackneg->GetSign() < 0) continue;
10eaad41	149
be1a7181	150	AliExternalTrackParam tpcTrackneg = (AliExternalTrackParam )trackneg->GetTPCInnerParam();
	151	if(!tpcTrackneg) continue;
	152	ptrackneg[0]=tpcTrackneg->Px();
	153	ptrackneg[1]=tpcTrackneg->Py();
	154	ptrackneg[2]=tpcTrackneg->Pz();
	155
	156	Bool_t firstLevelAcceptNegTrack=IsAcceptedForKink(esd, vertex, trackneg);
	157	if(firstLevelAcceptNegTrack==kFALSE) continue;
10eaad41	158
be1a7181	159	TVector3 negTrackMom(ptrackneg[0],ptrackneg[1],ptrackneg[2]);
10eaad41	160
be1a7181	161	Bool_t secondLevelAcceptNegTrack=IsAcceptedForTrack(esd, vertex, trackneg);
	162	if(secondLevelAcceptNegTrack==kFALSE) continue;
	163
	164	p4neg.SetVectM(negTrackMom, daughter2Mass);
10eaad41	165
	166	p4comb=p4pos;
	167	p4comb+=p4neg;
be1a7181	168
894840ad	169	if(p4comb.Vect().Pt()<=fMinPtTrackCut) continue;
be1a7181	170
894840ad	171	if((TMath::Abs(p4pos.Vect().Eta())<fuppereta)&&(TMath::Abs(p4neg.Vect().Eta())<fuppereta)&&(p4comb.Vect().Eta()<fuppereta)) {
be1a7181	172
	173	fPosKaonLikeSign->Fill(p4comb.M());
	174	fLikeSignInvmassPt->Fill(p4comb.M(), p4comb.Vect().Pt());
	175
	176	}
10eaad41	177
	178	} //inner track loop
	179
	180	} //outer track loop
	181
	182	PostData(1, fListOfHistos);
	183	}
	184
be1a7181	185	//________________________________________________________________________
	186	void AliResonanceKinkLikeSign::Terminate(Option_t *)
	187	{
	188
	189	}
	190
10eaad41	191	//____________________________________________________________________//
	192
	193	Float_t AliResonanceKinkLikeSign::GetSigmaToVertex(AliESDtrack* esdTrack) const {
	194	// Calculates the number of sigma to the vertex.
	195
	196	Float_t b[2];
	197	Float_t bRes[2];
	198	Float_t bCov[3];
	199
	200	esdTrack->GetImpactParameters(b,bCov);
	201
	202	if (bCov[0]<=0 \|\| bCov[2]<=0) {
	203	//AliDebug(1, "Estimated b resolution lower or equal zero!");
	204	bCov[0]=0; bCov[2]=0;
	205	}
	206	bRes[0] = TMath::Sqrt(bCov[0]);
	207	bRes[1] = TMath::Sqrt(bCov[2]);
	208
	209	if (bRes[0] == 0 \|\| bRes[1] ==0) return -1;
	210
	211	Float_t d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
	212
	213	if (TMath::Exp(-d * d / 2) < 1e-10) return 1000;
	214
	215	d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
	216
	217	return d;
	218	}
	219
	220	//____________________________________________________________________//
	221
	222	const AliESDVertex* AliResonanceKinkLikeSign::GetEventVertex(const AliESDEvent* esd) const
	223
	224	{
	225	// Get the vertex
	226
	227	const AliESDVertex* vertex = esd->GetPrimaryVertex();
	228
	229	if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
	230	else
	231	{
	232	vertex = esd->GetPrimaryVertexSPD();
	233	if((vertex->GetStatus()==kTRUE)&&(vertex->GetNContributors()>2)) return vertex;
	234	else
	235	return 0;
	236	}
	237	}
	238
be1a7181	239	//________________________________________________________________________
	240
	241	Bool_t AliResonanceKinkLikeSign::IsAcceptedForKink(AliESDEvent *localesd,
	242	const AliESDVertex localvertex, AliESDtrack localtrack) {
	243	// Apply the selections for each kink
	244
	245	Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
	246	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
	247	Double_t dca3D = 0.0;
	248
	249	AliExternalTrackParam tpcTrack = (AliExternalTrackParam )localtrack->GetTPCInnerParam();
	250	if(!tpcTrack) {
	251	gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
	252	dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
	253	cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
	254	}
	255	else {
	256	gPt = tpcTrack->Pt();
	257	gPx = tpcTrack->Px();
	258	gPy = tpcTrack->Py();
	259	gPz = tpcTrack->Pz();
	260	tpcTrack->PropagateToDCA(localvertex,
	261	localesd->GetMagneticField(),100.,dca,cov);
	262	}
	263
	264	if(GetSigmaToVertex(localtrack) > fMaxNSigmaToVertex) {
	265	if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a %lf sigmas to vertex TPC (max. requested: %lf)", GetSigmaToVertex(localtrack),fMaxNSigmaToVertex);
	266	return kFALSE;
	267	}
	268
	269	if(TMath::Abs(dca[0]) > fMaxDCAxy) {
	270	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);
	271	return kFALSE;
	272	}
	273
	274	if(TMath::Abs(dca[1]) > fMaxDCAzaxis) {
	275	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);
	276	return kFALSE;
	277	}
	278
	279	if(gPt <= fMinPtTrackCut) {
	280	if (fDebug > 1) Printf("IsAcceptedKink: Track rejected because it has a min value of pt of %lf (min. requested: %lf)", gPt, fMinPtTrackCut);
	281	return kFALSE;
	282	}
	283
	284	return kTRUE;
	285	}
	286
	287	//________________________________________________________________________
	288	Bool_t AliResonanceKinkLikeSign::IsAcceptedForTrack(AliESDEvent localesd, const AliESDVertex localvertex, AliESDtrack *localtrack) {
	289	// Apply the selections for each track
	290
	291	Double_t gPt = 0.0, gPx = 0.0, gPy = 0.0, gPz = 0.0;
	292	Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0}; //The impact parameters and their covariance.
	293	Double_t dca3D = 0.0;
	294
	295	AliExternalTrackParam tpcTrack = (AliExternalTrackParam )localtrack->GetTPCInnerParam();
	296	if(!tpcTrack) {
	297	gPt = 0.0; gPx = 0.0; gPy = 0.0; gPz = 0.0;
	298	dca[0] = -100.; dca[1] = -100.; dca3D = -100.;
	299	cov[0] = -100.; cov[1] = -100.; cov[2] = -100.;
	300	}
	301	else {
	302	gPt = tpcTrack->Pt();
303	gPx = tpcTrack->Px();
304	gPy = tpcTrack->Py();
305	gPz = tpcTrack->Pz();
306	tpcTrack->PropagateToDCA(localvertex,
307	localesd->GetMagneticField(),100.,dca,cov);
308	}
309
310	Int_t fcls[200];
311	Int_t nClustersTPC=localtrack->GetTPCclusters(fcls);
312	Float_t chi2perTPCcluster=-1.0;
313	if(nClustersTPC!=0) chi2perTPCcluster=(localtrack->GetTPCchi2())/Float_t(nClustersTPC);
314
315	Double_t extCov[15];
316	localtrack->GetExternalCovariance(extCov);
317
318	if((localtrack->GetStatus() & AliESDtrack::kTPCrefit) == 0) {
319	if (fDebug > 1) Printf("IsAccepted: Track rejected because of no refited in TPC");
320	return kFALSE;
321	}
322
323	if(nClustersTPC < fMinTPCclusters) {
324	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of nclusters (TPC) of %ld (min. requested: %ld)", nClustersTPC, fMinTPCclusters);
325	return kFALSE;
326	}
327
328	if(chi2perTPCcluster > fMaxChi2PerTPCcluster) {
329	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of chi2perTPCcluster of %lf (max. requested: %lf)", chi2perTPCcluster, fMaxChi2PerTPCcluster);
330	return kFALSE;
331	}
332
333	if(extCov[0] > fMaxCov0) {
334	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[0] of %lf (max. requested: %lf)", cov[0], fMaxCov0);
335	return kFALSE;
336	}
337
338	if(extCov[2] > fMaxCov2) {
339	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[2] of %lf (max. requested: %lf)", cov[2], fMaxCov2);
340	return kFALSE;
341	}
342
343	if(extCov[5] > fMaxCov5) {
344	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[5] of %lf (max. requested: %lf)", cov[5], fMaxCov5);
345	return kFALSE;
346	}
347
348	if(extCov[9] > fMaxCov9) {
349	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[9] of %lf (max. requested: %lf)", cov[9], fMaxCov9);
350	return kFALSE;
351	}
352
353	if(extCov[14] > fMaxCov14) {
354	if (fDebug > 1) Printf("IsAccepted: Track rejected because it has a value of cov[14] of %lf (max. requested: %lf)", cov[14], fMaxCov14);
355	return kFALSE;
356	}
357
358	return kTRUE;
359
360	}
361
362	//________________________________________________________________________
363	Bool_t AliResonanceKinkLikeSign::IsKink(AliESDEvent *localesd, Int_t kinkIndex, TVector3 trackMom)
364	{
365	// Test some kinematical criteria for each kink
366
367	AliESDkink *kink=localesd->GetKink(TMath::Abs(kinkIndex)-1);
368	const TVector3 motherMfromKink(kink->GetMotherP());
369	const TVector3 daughterMKink(kink->GetDaughterP());
370	Float_t qt=kink->GetQt();
371
372	Double_t maxDecAngKmu=f1->Eval(motherMfromKink.Mag(),0.,0.,0.);
373	Double_t maxDecAngpimu=f2->Eval(motherMfromKink.Mag(),0.,0.,0.);
374
375	Float_t kinkAngle=TMath::RadToDeg()*kink->GetAngle(2);
376
377	Float_t energyDaughterMu=TMath::Sqrt(daughterMKink.Mag()daughterMKink.Mag()+0.1056580.105658);
378	Float_t p1XM= motherMfromKink.Px();
379	Float_t p1YM= motherMfromKink.Py();
380	Float_t p1ZM= motherMfromKink.Pz();
381	Float_t p2XM= daughterMKink.Px();
382	Float_t p2YM= daughterMKink.Py();
383	Float_t p2ZM= daughterMKink.Pz();
384	Float_t p3Daughter=TMath::Sqrt(((p1XM-p2XM)(p1XM-p2XM))+((p1YM-p2YM)(p1YM-p2YM))+((p1ZM-p2ZM)*(p1ZM-p2ZM)));
385	Double_t invariantMassKmu= TMath::Sqrt((energyDaughterMu+p3Daughter)(energyDaughterMu+p3Daughter)-motherMfromKink.Mag()motherMfromKink.Mag());
386
894840ad	387	if((kinkAngle>maxDecAngpimu)&&(qt>fminQt)&&(qt<fmaxQt)&&((kink->GetR()>fminKinkRadius)&&(kink->GetR()<fmaxKinkRadius))&&(TMath::Abs(trackMom.Eta())<fuppereta)&&(invariantMassKmu<0.6)) {
be1a7181	388
	389	if(trackMom.Mag()<=1.1) {
	390	return kTRUE;
	391	}
	392	else
	393	if (kinkAngle<maxDecAngKmu) {
	394	return kTRUE;
	395	}
	396	}
	397	return kFALSE;
	398	}