[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsDstoKKpi.cxx

/**************************************************************************
 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * 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 for cuts on AOD reconstructed Ds->KKpi
//
// Author: A.Dainese, andrea.dainese@pd.infn.it
/////////////////////////////////////////////////////////////

#include <TDatabasePDG.h>
#include <Riostream.h>

#include "AliRDHFCutsDstoKKpi.h"
#include "AliAODRecoDecayHF3Prong.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"

ClassImp(AliRDHFCutsDstoKKpi)

//--------------------------------------------------------------------------
AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) : 
AliRDHFCuts(name)
{
  //
  // Default Constructor
  //
  Int_t nvars=14;
  SetNVars(nvars);
  TString varNames[14]={"inv. mass [GeV]",   
			"pTK [GeV/c]",
			"pTPi [GeV/c]",
			"d0K [cm]",
			"d0Pi [cm]",
			"dist12 [cm]",
			"sigmavert [cm]",
			"decLen [cm]",
			"ptMax [GeV/c]",
			"cosThetaPoint",
			"Sum d0^2 (cm^2)",
			"dca [cm]",
			"inv. mass (Mphi-MKK) [GeV]",
			"inv. mass (MKo*-MKpi) [GeV]"};
  Bool_t isUpperCut[14]={kTRUE,
			 kFALSE,
			 kFALSE,
			 kFALSE,
			 kFALSE,
			 kFALSE,
			 kTRUE,
			 kFALSE,
			 kFALSE,
			 kFALSE,
			 kFALSE,
			 kTRUE,
			 kTRUE,
			 kTRUE};
  SetVarNames(14,varNames,isUpperCut);
  Bool_t forOpt[14]={kFALSE,
		    kFALSE,
		    kFALSE,
		    kFALSE,
		    kFALSE,
		    kFALSE,
		    kTRUE,
		    kTRUE,
		    kTRUE,
		    kTRUE,
		    kTRUE,
		    kFALSE,
		    kTRUE,
		    kTRUE};
  SetVarsForOpt(7,forOpt);
  Float_t limits[2]={0,999999999.};
  SetPtBins(2,limits);
  if(fPidHF)delete fPidHF;
  fPidHF=new AliAODPidHF();
  Double_t plim[2]={0.6,0.8};
  Double_t nsigma[5]={2.,1.,2.,3.,0.};
  
  fPidHF->SetPLimit(plim);
  fPidHF->SetAsym(kTRUE);
  fPidHF->SetSigma(nsigma);
  fPidHF->SetMatch(1);
  fPidHF->SetTPC(1);
  fPidHF->SetTOF(1);
  fPidHF->SetITS(0);
  fPidHF->SetTRD(0);
  fPidHF->SetCompat(kTRUE);

}
//--------------------------------------------------------------------------
AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
  AliRDHFCuts(source)
{
  //
  // Copy constructor
  //

}
//--------------------------------------------------------------------------
AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
{
  //
  // assignment operator
  //
  if(&source == this) return *this;

  AliRDHFCuts::operator=(source);

  return *this;
}


//---------------------------------------------------------------------------
void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
  // 
  // Fills in vars the values of the variables 
  //

  if(nvars!=fnVarsForOpt) {
    printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
    return;
  }

  AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 
  Int_t iter=-1;
  if(fVarsForOpt[0]){
    iter++;
    if(TMath::Abs(pdgdaughters[0]==321)){
      vars[iter]=dd->InvMassDsKKpi();
    }else{
      vars[iter]=dd->InvMassDspiKK();
    }
  }
  if(fVarsForOpt[1]){
    iter++;
    Float_t minPtDau=99999.;
    for(Int_t iprong=0;iprong<3;iprong++){
      if(TMath::Abs(pdgdaughters[iprong])==321 && 
	 dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
    }
    vars[iter]=minPtDau;
  }
  if(fVarsForOpt[2]){
    iter++;
    for(Int_t iprong=0;iprong<3;iprong++){
      if(TMath::Abs(pdgdaughters[iprong])==211) {
	vars[iter]=dd->PtProng(iprong);
      }
    }
  }
  if(fVarsForOpt[3]){
    iter++;
    Float_t minImpParDau=99999.;
    for(Int_t iprong=0;iprong<3;iprong++){
      if(TMath::Abs(pdgdaughters[iprong])==321 &&
	 dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
    }
    vars[iter]=minImpParDau;
  }
  if(fVarsForOpt[4]){
    iter++;
    for(Int_t iprong=0;iprong<3;iprong++){
      if(TMath::Abs(pdgdaughters[iprong])==211) {
	vars[iter]=dd->Getd0Prong(iprong);
      }
    }
  }
  if(fVarsForOpt[5]){
    iter++;
    Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
    vars[iter]=minDistPair;
  }
  if(fVarsForOpt[6]){
    iter++;
    vars[iter]=dd->GetSigmaVert();
  }
  if(fVarsForOpt[7]){
    iter++;
    vars[iter] = dd->DecayLength();
  }
  if(fVarsForOpt[8]){
    iter++;
    Float_t ptmax=0;
    for(Int_t i=0;i<3;i++){
      if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
    }
    vars[iter]=ptmax;
  }
  if(fVarsForOpt[9]){
    iter++;
    vars[iter]=dd->CosPointingAngle();
  }
  if(fVarsForOpt[10]){
    iter++;
    vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
  }
  if(fVarsForOpt[11]){
    iter++;
    Float_t maxDCA=0.;
    for(Int_t i=0;i<3;i++){ 
      if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
    }
    vars[iter]=maxDCA;
  }
  if(fVarsForOpt[12]){
    iter++;
    Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
    if(TMath::Abs(pdgdaughters[0]==321)){
      
      Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
       vars[iter]=TMath::Abs(phimass01-mPDGPhi);
       // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
    }else{
      Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
       vars[iter]=TMath::Abs(phimass12-mPDGPhi);
       // vars[iter]=dd->InvMass2Prongs(1,2,321,321);      
    }
  }
  if(fVarsForOpt[13]){
    iter++;
    Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
    if(TMath::Abs(pdgdaughters[0]==321)){
      
      Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
      vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
      //	      vars[iter]=dd->InvMass2Prongs(1,2,321,211);
    }else{
      Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
      vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
      //	vars[iter]=dd->InvMass2Prongs(0,1,211,321);      
    }
  }

  
  return;
}
//---------------------------------------------------------------------------
Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
{
  //
  // Checking if Ds is in fiducial acceptance region 
  //

  if(pt > 5.) {
    // applying cut for pt > 5 GeV
    AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt)); 
    if (TMath::Abs(y) > 0.8) return kFALSE;
    
  } else {
    // appliying smooth cut for pt < 5 GeV
    Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; 
    Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;		
    AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); 
    if (y < minFiducialY || y > maxFiducialY) return kFALSE;    
  }

  return kTRUE;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
  // PID selection
  // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both 
  Int_t retCode=3;
  Bool_t okKKpi=kTRUE;
  Bool_t okpiKK=kTRUE;
  if(!fUsePID || !rd) return retCode;
  if(!fPidHF){
    AliWarning("AliAODPidHF not created!");
    return retCode;
  }
  Int_t nKaons=0;
  Int_t nNotKaons=0;
  Int_t sign= rd->GetCharge(); 
  for(Int_t iDaught=0; iDaught<3; iDaught++){
    AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
    Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
    Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
    Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
    
    if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
    if(sign!=track->Charge()){// must be kaon
      if(isKaon<0) return 0;
    }
    if(isKaon>0 && isPion<0) nKaons++;
    if(isKaon<0) nNotKaons++;
    if(iDaught==0){
      if(isKaon<0) okKKpi=kFALSE;
      if(isPion<0) okpiKK=kFALSE;
    }
    else if(iDaught==2){
      if(isKaon<0) okpiKK=kFALSE;
      if(isPion<0) okKKpi=kFALSE;
    }
  }
  
  if(nKaons>2)return 0;
  if(nNotKaons>1) return 0;
  
  if(!okKKpi) retCode-=1;
  if(!okpiKK) retCode-=2;

  return retCode;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
  //
  // Apply selection
  //

  if(!fCutsRD){
    cout<<"Cut matrice not inizialized. Exit..."<<endl;
    return 0;
  }
  //PrintAll();
  AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;

  if(!d){
    cout<<"AliAODRecoDecayHF3Prong null"<<endl;
    return 0;
  }


  // selection on daughter tracks 
  if(selectionLevel==AliRDHFCuts::kAll || 
     selectionLevel==AliRDHFCuts::kTracks) {
    if(!AreDaughtersSelected(d)) return 0;
  }


  // PID selection
  Int_t returnvaluePID=3;  
  if(selectionLevel==AliRDHFCuts::kAll || 
     selectionLevel==AliRDHFCuts::kCandidate ||     
     selectionLevel==AliRDHFCuts::kPID) {
    returnvaluePID = IsSelectedPID(d);
  }
  if(returnvaluePID==0)return 0;
  Bool_t okPidDsKKpi=returnvaluePID&1;
  Bool_t okPidDspiKK=returnvaluePID&2;
 
  Int_t returnvalueCuts=15;
  // selection on candidate
  if(selectionLevel==AliRDHFCuts::kAll || 
     selectionLevel==AliRDHFCuts::kCandidate) {
    //recalculate vertex w/o daughters
    AliAODVertex *origownvtx=0x0;
    AliAODVertex *recvtx=0x0;
    if(fRemoveDaughtersFromPrimary) {
      if(!aod) {
	AliError("Can not remove daughters from vertex without AOD event");
	return 0;
      }
      if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
      recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
      if(!recvtx){
	AliDebug(2,"Removal of daughter tracks failed");
	//recvtx=d->GetPrimaryVtx();
	if(origownvtx){
	  delete origownvtx;
	  origownvtx=NULL;
	}
	return 0;
      }
      //set recalculed primary vertex
      d->SetOwnPrimaryVtx(recvtx);
      delete recvtx; recvtx=NULL;
    }

    Int_t okDsKKpi=1;
    Int_t okDspiKK=1;
    Int_t okMassPhi=0;
    Int_t okMassK0star=0;

    Double_t pt=d->Pt();
    Int_t ptbin=PtBin(pt);

    Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
    Double_t mDsKKpi=d->InvMassDsKKpi();
    Double_t mDspiKK=d->InvMassDspiKK();
    if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
    if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
    if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;
    if(okPidDsKKpi && !okDsKKpi)  returnvalueCuts=0;
    if(okPidDspiKK && !okDspiKK) returnvalueCuts=0;

    //single track
    if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || 
       TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) returnvalueCuts=0;
    if(okDsKKpi){
      if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(1,ptbin)] || 
	 TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
      if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] || 
	 TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
    }
    if(okDspiKK){
      if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] || 
	 TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
      if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(1,ptbin)] || 
	 TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
    }
    if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;

    // cuts on resonant decays (via Phi or K0*)
    Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
    Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
    if(okDsKKpi){
      Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
      Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
      if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
      if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
      if(!okMassPhi && !okMassK0star) okDsKKpi=0;
    }
    if(okDspiKK){
      Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
      Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
      if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
      if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
      if(!okMassPhi && !okMassK0star) okDspiKK=0;
    }
    if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;

    // Cuts on track pairs
    for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)])  returnvalueCuts=0;
    if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] || 
       d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) returnvalueCuts=0;


    // Cuts on candidate triplet
    if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]) returnvalueCuts=0;
    if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) returnvalueCuts=0;
    if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] && 
       TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] && 
       TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)]) returnvalueCuts=0;
    if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)])returnvalueCuts=0;
    Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
    if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)])returnvalueCuts=0;

     // unset recalculated primary vertex when not needed any more
    if(origownvtx) {
      d->SetOwnPrimaryVtx(origownvtx);
      delete origownvtx;
      origownvtx=NULL;
    } else if(fRemoveDaughtersFromPrimary) {
      d->UnsetOwnPrimaryVtx();
      AliDebug(3,"delete new vertex\n");
    }
   
    if(returnvalueCuts == 0) return 0;
    Int_t returnvalue=0;
    if(okDsKKpi) returnvalue+=1;
    if(okDspiKK) returnvalue+=2;
    if(okMassPhi) returnvalue+=4;
    if(okMassK0star) returnvalue+=8;

    return returnvalue;
  }
  return returnvalueCuts;

}
//---------------------------------------------------------------------------
Commit	Line	Data
e3d40058	1	/**************************************************************************
	2	* Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/////////////////////////////////////////////////////////////
	17	//
	18	// Class for cuts on AOD reconstructed Ds->KKpi
	19	//
	20	// Author: A.Dainese, andrea.dainese@pd.infn.it
	21	/////////////////////////////////////////////////////////////
	22
	23	#include <TDatabasePDG.h>
	24	#include <Riostream.h>
	25
	26	#include "AliRDHFCutsDstoKKpi.h"
	27	#include "AliAODRecoDecayHF3Prong.h"
	28	#include "AliAODTrack.h"
	29	#include "AliESDtrack.h"
	30
	31	ClassImp(AliRDHFCutsDstoKKpi)
	32
	33	//--------------------------------------------------------------------------
a9b75906	34	AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) :
a9b75906	35	AliRDHFCuts(name)
e3d40058	36	{
	37	//
	38	// Default Constructor
	39	//
a9b75906	40	Int_t nvars=14;
e3d40058	41	SetNVars(nvars);
a9b75906	42	TString varNames[14]={"inv. mass [GeV]",
4755453e	43	"pTK [GeV/c]",
	44	"pTPi [GeV/c]",
	45	"d0K [cm]",
	46	"d0Pi [cm]",
	47	"dist12 [cm]",
	48	"sigmavert [cm]",
	49	"decLen [cm]",
	50	"ptMax [GeV/c]",
	51	"cosThetaPoint",
	52	"Sum d0^2 (cm^2)",
	53	"dca [cm]",
	54	"inv. mass (Mphi-MKK) [GeV]",
a9b75906	55	"inv. mass (MKo*-MKpi) [GeV]"};
a9b75906	56	Bool_t isUpperCut[14]={kTRUE,
4755453e	57	kFALSE,
	58	kFALSE,
	59	kFALSE,
	60	kFALSE,
	61	kFALSE,
	62	kTRUE,
	63	kFALSE,
	64	kFALSE,
	65	kFALSE,
	66	kFALSE,
	67	kTRUE,
	68	kTRUE,
	69	kTRUE};
a9b75906	70	SetVarNames(14,varNames,isUpperCut);
a9b75906	71	Bool_t forOpt[14]={kFALSE,
4755453e	72	kFALSE,
e3d40058	73	kFALSE,
	74	kFALSE,
	75	kFALSE,
	76	kFALSE,
	77	kTRUE,
4755453e	78	kTRUE,
	79	kTRUE,
	80	kTRUE,
	81	kTRUE,
	82	kFALSE,
	83	kTRUE,
e3d40058	84	kTRUE};
4755453e	85	SetVarsForOpt(7,forOpt);
e3d40058	86	Float_t limits[2]={0,999999999.};
e3d40058	87	SetPtBins(2,limits);
3787d195	88	if(fPidHF)delete fPidHF;
	89	fPidHF=new AliAODPidHF();
	90	Double_t plim[2]={0.6,0.8};
	91	Double_t nsigma[5]={2.,1.,2.,3.,0.};
	92
	93	fPidHF->SetPLimit(plim);
	94	fPidHF->SetAsym(kTRUE);
	95	fPidHF->SetSigma(nsigma);
	96	fPidHF->SetMatch(1);
	97	fPidHF->SetTPC(1);
	98	fPidHF->SetTOF(1);
	99	fPidHF->SetITS(0);
	100	fPidHF->SetTRD(0);
	101	fPidHF->SetCompat(kTRUE);
	102
e3d40058	103	}
	104	//--------------------------------------------------------------------------
	105	AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
	106	AliRDHFCuts(source)
	107	{
	108	//
	109	// Copy constructor
	110	//
	111
	112	}
	113	//--------------------------------------------------------------------------
	114	AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
	115	{
	116	//
	117	// assignment operator
	118	//
	119	if(&source == this) return *this;
	120
	121	AliRDHFCuts::operator=(source);
	122
	123	return *this;
	124	}
	125
	126
	127	//---------------------------------------------------------------------------
	128	void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF d,Float_t vars,Int_t nvars,Int_t *pdgdaughters) {
	129	//
	130	// Fills in vars the values of the variables
	131	//
	132
	133	if(nvars!=fnVarsForOpt) {
	134	printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
	135	return;
	136	}
	137
	138	AliAODRecoDecayHF3Prong dd = (AliAODRecoDecayHF3Prong)d;
	139
4755453e	140	Int_t iter=-1;
	141	if(fVarsForOpt[0]){
	142	iter++;
	143	if(TMath::Abs(pdgdaughters[0]==321)){
	144	vars[iter]=dd->InvMassDsKKpi();
	145	}else{
	146	vars[iter]=dd->InvMassDspiKK();
	147	}
	148	}
	149	if(fVarsForOpt[1]){
	150	iter++;
	151	Float_t minPtDau=99999.;
	152	for(Int_t iprong=0;iprong<3;iprong++){
	153	if(TMath::Abs(pdgdaughters[iprong])==321 &&
	154	dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
	155	}
	156	vars[iter]=minPtDau;
	157	}
	158	if(fVarsForOpt[2]){
	159	iter++;
	160	for(Int_t iprong=0;iprong<3;iprong++){
	161	if(TMath::Abs(pdgdaughters[iprong])==211) {
	162	vars[iter]=dd->PtProng(iprong);
	163	}
	164	}
	165	}
	166	if(fVarsForOpt[3]){
	167	iter++;
	168	Float_t minImpParDau=99999.;
	169	for(Int_t iprong=0;iprong<3;iprong++){
	170	if(TMath::Abs(pdgdaughters[iprong])==321 &&
	171	dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
	172	}
	173	vars[iter]=minImpParDau;
	174	}
	175	if(fVarsForOpt[4]){
	176	iter++;
	177	for(Int_t iprong=0;iprong<3;iprong++){
	178	if(TMath::Abs(pdgdaughters[iprong])==211) {
	179	vars[iter]=dd->Getd0Prong(iprong);
	180	}
	181	}
	182	}
	183	if(fVarsForOpt[5]){
	184	iter++;
	185	Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
	186	vars[iter]=minDistPair;
	187	}
	188	if(fVarsForOpt[6]){
	189	iter++;
	190	vars[iter]=dd->GetSigmaVert();
	191	}
	192	if(fVarsForOpt[7]){
	193	iter++;
	194	vars[iter] = dd->DecayLength();
	195	}
	196	if(fVarsForOpt[8]){
	197	iter++;
	198	Float_t ptmax=0;
	199	for(Int_t i=0;i<3;i++){
	200	if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
	201	}
	202	vars[iter]=ptmax;
	203	}
204	if(fVarsForOpt[9]){
205	iter++;
206	vars[iter]=dd->CosPointingAngle();
207	}
208	if(fVarsForOpt[10]){
209	iter++;
210	vars[iter]=dd->Getd0Prong(0)dd->Getd0Prong(0)+dd->Getd0Prong(1)dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
211	}
212	if(fVarsForOpt[11]){
213	iter++;
214	Float_t maxDCA=0.;
215	for(Int_t i=0;i<3;i++){
216	if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
217	}
218	vars[iter]=maxDCA;
219	}
220	if(fVarsForOpt[12]){
221	iter++;
2c3decce	222	Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
4755453e	223	if(TMath::Abs(pdgdaughters[0]==321)){
2c3decce	224
	225	Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
	226	vars[iter]=TMath::Abs(phimass01-mPDGPhi);
	227	// vars[iter]=dd->InvMass2Prongs(0,1,321,321);
4755453e	228	}else{
2c3decce	229	Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
	230	vars[iter]=TMath::Abs(phimass12-mPDGPhi);
	231	// vars[iter]=dd->InvMass2Prongs(1,2,321,321);
4755453e	232	}
	233	}
	234	if(fVarsForOpt[13]){
	235	iter++;
2c3decce	236	Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
4755453e	237	if(TMath::Abs(pdgdaughters[0]==321)){
2c3decce	238
	239	Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
	240	vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
	241	// vars[iter]=dd->InvMass2Prongs(1,2,321,211);
4755453e	242	}else{
2c3decce	243	Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
	244	vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
	245	// vars[iter]=dd->InvMass2Prongs(0,1,211,321);
4755453e	246	}
	247	}
	248
e3d40058	249
	250	return;
	251	}
3787d195	252	//---------------------------------------------------------------------------
	253	Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
	254	{
	255	//
	256	// Checking if Ds is in fiducial acceptance region
	257	//
	258
	259	if(pt > 5.) {
	260	// applying cut for pt > 5 GeV
	261	AliDebug(2,Form("pt of Ds = %f (> 5), cutting at \|y\| < 0.8",pt));
	262	if (TMath::Abs(y) > 0.8) return kFALSE;
	263
	264	} else {
	265	// appliying smooth cut for pt < 5 GeV
	266	Double_t maxFiducialY = -0.2/15ptpt+1.9/15*pt+0.5;
	267	Double_t minFiducialY = 0.2/15ptpt-1.9/15*pt-0.5;
	268	AliDebug(2,Form("pt of Ds = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
	269	if (y < minFiducialY \|\| y > maxFiducialY) return kFALSE;
	270	}
	271
	272	return kTRUE;
	273	}
	274
45ab4739	275	//---------------------------------------------------------------------------
	276	Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
	277	// PID selection
	278	// return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both
	279	Int_t retCode=3;
	280	Bool_t okKKpi=kTRUE;
	281	Bool_t okpiKK=kTRUE;
	282	if(!fUsePID \|\| !rd) return retCode;
	283	if(!fPidHF){
	284	AliWarning("AliAODPidHF not created!");
	285	return retCode;
	286	}
	287	Int_t nKaons=0;
	288	Int_t nNotKaons=0;
	289	Int_t sign= rd->GetCharge();
	290	for(Int_t iDaught=0; iDaught<3; iDaught++){
	291	AliAODTrack track=(AliAODTrack)rd->GetDaughter(iDaught);
	292	Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
	293	Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
	294	Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
	295
	296	if(isProton>0 && isKaon<0 && isPion<0) return 0;
	297	if(sign!=track->Charge()){// must be kaon
	298	if(isKaon<0) return 0;
	299	}
	300	if(isKaon>0 && isPion<0) nKaons++;
	301	if(isKaon<0) nNotKaons++;
	302	if(iDaught==0){
	303	if(isKaon<0) okKKpi=kFALSE;
	304	if(isPion<0) okpiKK=kFALSE;
	305	}
	306	else if(iDaught==2){
	307	if(isKaon<0) okpiKK=kFALSE;
	308	if(isPion<0) okKKpi=kFALSE;
	309	}
	310	}
	311
	312	if(nKaons>2)return 0;
	313	if(nNotKaons>1) return 0;
	314
	315	if(!okKKpi) retCode-=1;
	316	if(!okpiKK) retCode-=2;
	317
	318	return retCode;
	319	}
	320
e3d40058	321	//---------------------------------------------------------------------------
3787d195	322	Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
e3d40058	323	//
	324	// Apply selection
	325	//
	326
	327	if(!fCutsRD){
	328	cout<<"Cut matrice not inizialized. Exit..."<<endl;
	329	return 0;
	330	}
	331	//PrintAll();
	332	AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
	333
	334	if(!d){
	335	cout<<"AliAODRecoDecayHF3Prong null"<<endl;
	336	return 0;
	337	}
	338
	339
	340	// selection on daughter tracks
	341	if(selectionLevel==AliRDHFCuts::kAll \|\|
	342	selectionLevel==AliRDHFCuts::kTracks) {
	343	if(!AreDaughtersSelected(d)) return 0;
	344	}
	345
	346
e3d40058	347
45ab4739	348	// PID selection
	349	Int_t returnvaluePID=3;
	350	if(selectionLevel==AliRDHFCuts::kAll \|\|
	351	selectionLevel==AliRDHFCuts::kCandidate \|\|
	352	selectionLevel==AliRDHFCuts::kPID) {
	353	returnvaluePID = IsSelectedPID(d);
	354	}
	355	if(returnvaluePID==0)return 0;
	356	Bool_t okPidDsKKpi=returnvaluePID&1;
	357	Bool_t okPidDspiKK=returnvaluePID&2;
	358
3787d195	359	Int_t returnvalueCuts=15;
e3d40058	360	// selection on candidate
	361	if(selectionLevel==AliRDHFCuts::kAll \|\|
	362	selectionLevel==AliRDHFCuts::kCandidate) {
3787d195	363	//recalculate vertex w/o daughters
	364	AliAODVertex *origownvtx=0x0;
	365	AliAODVertex *recvtx=0x0;
	366	if(fRemoveDaughtersFromPrimary) {
	367	if(!aod) {
	368	AliError("Can not remove daughters from vertex without AOD event");
	369	return 0;
	370	}
	371	if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
	372	recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
	373	if(!recvtx){
	374	AliDebug(2,"Removal of daughter tracks failed");
	375	//recvtx=d->GetPrimaryVtx();
	376	if(origownvtx){
	377	delete origownvtx;
	378	origownvtx=NULL;
	379	}
	380	return 0;
	381	}
	382	//set recalculed primary vertex
	383	d->SetOwnPrimaryVtx(recvtx);
	384	delete recvtx; recvtx=NULL;
	385	}
e3d40058	386
4755453e	387	Int_t okDsKKpi=1;
	388	Int_t okDspiKK=1;
	389	Int_t okMassPhi=0;
	390	Int_t okMassK0star=0;
e3d40058	391
4755453e	392	Double_t pt=d->Pt();
	393	Int_t ptbin=PtBin(pt);
	394
	395	Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
	396	Double_t mDsKKpi=d->InvMassDsKKpi();
	397	Double_t mDspiKK=d->InvMassDspiKK();
	398	if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
	399	if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
3787d195	400	if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;
	401	if(okPidDsKKpi && !okDsKKpi) returnvalueCuts=0;
	402	if(okPidDspiKK && !okDspiKK) returnvalueCuts=0;
4755453e	403
	404	//single track
	405	if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] \|\|
3787d195	406	TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) returnvalueCuts=0;
4755453e	407	if(okDsKKpi){
	408	if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(1,ptbin)] \|\|
	409	TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
	410	if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] \|\|
	411	TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
	412	}
	413	if(okDspiKK){
	414	if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] \|\|
	415	TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
	416	if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(1,ptbin)] \|\|
	417	TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
	418	}
3787d195	419	if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;
4755453e	420
	421	// cuts on resonant decays (via Phi or K0*)
	422	Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
	423	Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
	424	if(okDsKKpi){
	425	Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
	426	Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
	427	if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
	428	if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
	429	if(!okMassPhi && !okMassK0star) okDsKKpi=0;
	430	}
	431	if(okDspiKK){
	432	Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
	433	Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
	434	if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
	435	if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
	436	if(!okMassPhi && !okMassK0star) okDspiKK=0;
	437	}
3787d195	438	if(!okDsKKpi && !okDspiKK) returnvalueCuts=0;
4755453e	439
4755453e	440	// Cuts on track pairs
3787d195	441	for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) returnvalueCuts=0;
4755453e	442	if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] \|\|
3787d195	443	d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) returnvalueCuts=0;
4755453e	444
	445
	446	// Cuts on candidate triplet
3787d195	447	if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]) returnvalueCuts=0;
3787d195	448	if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) returnvalueCuts=0;
4755453e	449	if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] &&
4755453e	450	TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] &&
3787d195	451	TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)]) returnvalueCuts=0;
3787d195	452	if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)])returnvalueCuts=0;
4755453e	453	Double_t sum2=d->Getd0Prong(0)d->Getd0Prong(0)+d->Getd0Prong(1)d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
3787d195	454	if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)])returnvalueCuts=0;
	455
	456	// unset recalculated primary vertex when not needed any more
	457	if(origownvtx) {
	458	d->SetOwnPrimaryVtx(origownvtx);
	459	delete origownvtx;
	460	origownvtx=NULL;
	461	} else if(fRemoveDaughtersFromPrimary) {
	462	d->UnsetOwnPrimaryVtx();
	463	AliDebug(3,"delete new vertex\n");
	464	}
	465
	466	if(returnvalueCuts == 0) return 0;
	467	Int_t returnvalue=0;
4755453e	468	if(okDsKKpi) returnvalue+=1;
	469	if(okDspiKK) returnvalue+=2;
	470	if(okMassPhi) returnvalue+=4;
	471	if(okMassK0star) returnvalue+=8;
	472
3787d195	473	return returnvalue;
4755453e	474	}
3787d195	475	return returnvalueCuts;
e3d40058	476
	477	}
	478	//---------------------------------------------------------------------------