[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliRDHFCutsLctopK0sfromAODtracks.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.                  *
 **************************************************************************/

/* $Id$ */

/////////////////////////////////////////////////////////////
//
// Class for cuts on AOD reconstructed Lc->p+K0s
//
// Modified by Y.S Watanabe - wyosuke@cns.s.u-tokyo.ac.jp
//
/////////////////////////////////////////////////////////////

#include <Riostream.h>

#include <TDatabasePDG.h>
#include <TMath.h>

#include "AliAnalysisManager.h"
#include "AliInputEventHandler.h"
#include "AliPIDResponse.h"
#include "AliRDHFCutsLctopK0sfromAODtracks.h"
#include "AliAODRecoCascadeHF.h"
#include "AliAODTrack.h"
#include "AliESDtrack.h"
#include "AliESDVertex.h"
#include "AliAODVertex.h"
#include "AliAODv0.h"
#include "AliAODcascade.h"
#include "AliESDv0.h"

using std::cout;
using std::endl;

ClassImp(AliRDHFCutsLctopK0sfromAODtracks)

//--------------------------------------------------------------------------
AliRDHFCutsLctopK0sfromAODtracks::AliRDHFCutsLctopK0sfromAODtracks(const char* name) :
AliRDHFCuts(name),
  fPIDStrategy(kNSigmaCuts),
  fCombinedPIDThreshold(0.),
  fUseOnTheFlyV0(kFALSE),
  fProdTrackPtMin(0.3),
  fProdTrackEtaRange(0.8),
  fProdUseAODFilterBit(kTRUE),
  fProdV0MassTolK0s(0.01),
  fProdV0PtMin(0.5),
  fProdV0CosPointingAngleToPrimVtxMin(0.99),
  fProdV0DcaDaughtersMax(1.5),
  fProdV0DaughterEtaRange(0.8),
  fProdV0DaughterPtMin(0.0),
  fProdV0DaughterTPCClusterMin(70),
fProdRoughMassTol(0.25),
fProdRoughPtMin(0.0),
fnCuts(0),
  fnTotalCutBins(1),
  fCutsArray(0)
{
  //
  // Default Constructor
  //

  const Int_t nvars=7;
  SetNVars(nvars);
  TString varNames[nvars]={"Lc inv. mass [GeV/c2]",                   //  0
			   "Lc pT [GeV/c]", //1
			   "Bachelor pT [GeV/c]", //2
			   "Bachelor d0 [cm]", //3
			   "V0 d0 [cm]", //4
			   "K0s mass [GeV/c2]", //5
			   "Decay Length XY [cm]" //6
  };

  Bool_t isUpperCut[nvars]={kTRUE,  //  0
			    kFALSE, //1
			    kFALSE, //2
			    kTRUE, //3
			    kTRUE, //4
			    kTRUE, //5
			    kFALSE //6
  };
  SetVarNames(nvars,varNames,isUpperCut);
  Bool_t forOpt[nvars]={kFALSE, //  0
			kFALSE, //1
			kTRUE, //2
			kTRUE, //3
			kTRUE, //4
			kTRUE, //5
			kTRUE //6
  };
  SetVarsForOpt(nvars,forOpt);

  Float_t limits[2]={0,999999999.};
  SetPtBins(2,limits);
}
//--------------------------------------------------------------------------
AliRDHFCutsLctopK0sfromAODtracks::AliRDHFCutsLctopK0sfromAODtracks(const AliRDHFCutsLctopK0sfromAODtracks &source) :
  AliRDHFCuts(source),
  fPIDStrategy(source.fPIDStrategy),
  fCombinedPIDThreshold(source.fCombinedPIDThreshold),
  fUseOnTheFlyV0(source.fUseOnTheFlyV0),
  fProdTrackPtMin(source.fProdTrackPtMin),
  fProdTrackEtaRange(source.fProdTrackEtaRange),
  fProdUseAODFilterBit(source.fProdUseAODFilterBit),
  fProdV0MassTolK0s(source.fProdV0MassTolK0s),
  fProdV0PtMin(source.fProdV0PtMin),
  fProdV0CosPointingAngleToPrimVtxMin(source.fProdV0CosPointingAngleToPrimVtxMin),
  fProdV0DcaDaughtersMax(source.fProdV0DcaDaughtersMax),
  fProdV0DaughterEtaRange(source.fProdV0DaughterEtaRange),
  fProdV0DaughterPtMin(source.fProdV0DaughterPtMin),
  fProdV0DaughterTPCClusterMin(source.fProdV0DaughterTPCClusterMin),
  fProdRoughMassTol(source.fProdRoughMassTol),
  fProdRoughPtMin(source.fProdRoughPtMin),
  fnCuts(source.fnCuts),
  fnTotalCutBins(source.fnTotalCutBins),
  fCutsArray(0)
{
  //
  // Copy constructor
  //
  if(source.fCutsArray) SetCutsArray(source.fnTotalCutBins,source.fCutsArray);
}
//--------------------------------------------------------------------------
AliRDHFCutsLctopK0sfromAODtracks &AliRDHFCutsLctopK0sfromAODtracks::operator=(const AliRDHFCutsLctopK0sfromAODtracks &source)
{
  //
  // assignment operator
  //

  if (this != &source) {
    AliRDHFCuts::operator=(source);
  }

  fPIDStrategy = source.fPIDStrategy;
  fCombinedPIDThreshold = source.fCombinedPIDThreshold;
  fUseOnTheFlyV0 = source.fUseOnTheFlyV0;
  fProdTrackPtMin = source.fProdTrackPtMin;
  fProdTrackEtaRange = source.fProdTrackEtaRange;
  fProdUseAODFilterBit = source.fProdUseAODFilterBit;
  fProdV0MassTolK0s = source.fProdV0MassTolK0s;
  fProdV0PtMin = source.fProdV0PtMin;
  fProdV0CosPointingAngleToPrimVtxMin = source.fProdV0CosPointingAngleToPrimVtxMin;
  fProdV0DcaDaughtersMax=source.fProdV0DcaDaughtersMax;
  fProdV0DaughterEtaRange=source.fProdV0DaughterEtaRange;
  fProdV0DaughterPtMin=source.fProdV0DaughterPtMin;
  fProdV0DaughterTPCClusterMin=source.fProdV0DaughterTPCClusterMin;
  fProdRoughMassTol = source.fProdRoughMassTol;
  fProdRoughPtMin = source.fProdRoughPtMin;
  
  if(source.fnCuts) fnCuts = source.fnCuts;
  if(source.fnTotalCutBins) fnTotalCutBins = source.fnCuts;
  if(source.fCutsArray) SetCutsArray(source.fnTotalCutBins,source.fCutsArray);
  
  return *this;
}

//---------------------------------------------------------------------------
AliRDHFCutsLctopK0sfromAODtracks::~AliRDHFCutsLctopK0sfromAODtracks() {
  //
  //  Default Destructor
  //
  
  if(fCutsArray) delete[] fCutsArray;
}

//---------------------------------------------------------------------------
void AliRDHFCutsLctopK0sfromAODtracks::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
  //
  // Fills in vars the values of the variables
  //
  
  if (pdgdaughters[0]==-9999) return; // dummy

  AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)d;
  if(!dd){
    AliError("No AliAODRecoCascadeHF object found\n");
    return;
  }
  
  if (nvars!=fnVarsForOpt) {
    AliError("Cut object seems to have the wrong number of variables\n");
    return;
  }

  //Double_t ptD=d->Pt();
  //Int_t ptbin=PtBin(ptD);
  Int_t iter=-1;

  if(fVarsForOpt[0]){
    iter++;
    Double_t mlcPDG =  TDatabasePDG::Instance()->GetParticle(4122)->Mass();
    vars[iter]= TMath::Abs(dd->InvMassLctoK0sP()-mlcPDG) ;
  }
  if(fVarsForOpt[1]){
    iter++;
    vars[iter]= dd->Pt();
  }
  if(fVarsForOpt[2]){
    iter++;
    AliAODTrack *part = dd->GetBachelor();
    vars[iter]= part->Pt();
  }
  if(fVarsForOpt[3]){
    iter++;
    vars[iter]= dd->Getd0Prong(0);
  }
  if(fVarsForOpt[4]){
    iter++;
    vars[iter]= dd->Getd0Prong(1);
  }
  if(fVarsForOpt[5]){
    iter++;
    AliAODv0 *v0 = dd->Getv0();
    vars[iter]= v0->MassK0Short();
  }
  if(fVarsForOpt[6]){
    iter++;
    vars[iter]= dd->DecayLengthXY();
  }
  
  return;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelected(TObject* obj, Int_t selectionLevel) {
  //
  // Apply selection
  //

  if (!fCutsRD) {
    AliFatal("Cut matrix not inizialized. Exit...");
    return 0;
  }

  AliAODRecoCascadeHF* d=(AliAODRecoCascadeHF*)obj;
  if(!d){
    AliDebug(2,"AliAODRecoCascadeHF null");
    return 0;
  }

  Double_t ptD=d->Pt();
  if(ptD<fMinPtCand) return 0;
  if(ptD>fMaxPtCand) return 0;

  if (selectionLevel==AliRDHFCuts::kAll ||
      selectionLevel==AliRDHFCuts::kTracks) {
    //Performed in production stage
  }

  Int_t returnvalueCuts=1;
  // selection on candidate
  if (selectionLevel==AliRDHFCuts::kAll ||
      selectionLevel==AliRDHFCuts::kCandidate) {
    
    Double_t pt=d->Pt();
    Int_t ptbin=PtBin(pt);
    if (ptbin==-1) {
      return 0;
    }
    Bool_t okcand=kTRUE;
    
    Double_t mlcPDG =  TDatabasePDG::Instance()->GetParticle(4122)->Mass();
    Double_t mk0sPDG =  TDatabasePDG::Instance()->GetParticle(310)->Mass();
    AliAODTrack *part = d->GetBachelor();
    AliAODv0 *v0 = d->Getv0();
    
    if(TMath::Abs(d->InvMassLctoK0sP()-mlcPDG) > fCutsRD[GetGlobalIndex(0,ptbin)])
      {
	okcand = kFALSE;
      }
    if(d->Pt() < fCutsRD[GetGlobalIndex(1,ptbin)])
      {
	okcand = kFALSE;
      }
    if(part->Pt() < fCutsRD[GetGlobalIndex(2,ptbin)])
      {
	okcand = kFALSE;
      }
    if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(3,ptbin)])
      {
	okcand = kFALSE;
      }
    if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(4,ptbin)])
      {
	okcand = kFALSE;
      }
    if(TMath::Abs(v0->MassK0Short()-mk0sPDG) > fCutsRD[GetGlobalIndex(5,ptbin)])
      {
	okcand = kFALSE;
      }
    if(d->DecayLengthXY() < fCutsRD[GetGlobalIndex(6,ptbin)])
      {
	okcand = kFALSE;
      }
    
    if(!okcand)  return 0;
    returnvalueCuts = 1;
  }
  
  Int_t returnvaluePID=1;
  if(selectionLevel==AliRDHFCuts::kAll ||
     selectionLevel==AliRDHFCuts::kCandidate|| 
     selectionLevel==AliRDHFCuts::kPID) {

    switch(fPIDStrategy){
    case kNSigmaCuts:
      returnvaluePID = IsSelectedPID(d);
      break;
    case kCombinedCuts:
      returnvaluePID = IsSelectedCombinedPID(d);
      break;
    }
  }
  
  Int_t returnvalue = 0;
  if(returnvalueCuts==1 && returnvaluePID==1) returnvalue=1;
  
  return returnvalue;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelectedPID(AliAODRecoDecayHF* obj) 
{
  //
  // IsSelectedPID
  //

  if(!fUsePID || !obj) return 1;

  AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)obj;
  AliAODTrack *part = dd->GetBachelor();
  
  Int_t returnvalue=1;

  if(fPidHF->GetPidResponse()==0x0){
    AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
    AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
    AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
    fPidHF->SetPidResponse(pidResp);
  }

  Int_t isProton=fPidHF->MakeRawPid(part,4); 
  if(isProton<1) returnvalue = 0;
  
  return returnvalue;
}

//---------------------------------------------------------------------------
Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
  //
  // IsSelectedCombinedPID
  //
    
  if(!fUsePID || !obj) {return 1;}

  AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)obj;
  AliAODTrack *part = dd->GetBachelor();
  if(!part) return 0;

  Int_t returnvalue=1;
  Double_t probProton = GetProtonProbabilityTPCTOF(part);
  if(probProton<fCombinedPIDThreshold) returnvalue = 0;

  return returnvalue;
}

//________________________________________________________________________
Double_t AliRDHFCutsLctopK0sfromAODtracks::GetProtonProbabilityTPCTOF(AliAODTrack* trk) 
{
  //
  // Get Proton probability
  //
  fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
  Double_t prob1[AliPID::kSPECIES];
  UInt_t detUsed1 = fPidHF->GetPidCombined()->ComputeProbabilities(trk, fPidHF->GetPidResponse(), prob1);
  if (detUsed1 != (UInt_t)fPidHF->GetPidCombined()->GetDetectorMask() ) 
    { 
      fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
      detUsed1 = fPidHF->GetPidCombined()->ComputeProbabilities(trk, fPidHF->GetPidResponse(), prob1);
    }
  return prob1[AliPID::kProton];
}

//________________________________________________________________________
Bool_t AliRDHFCutsLctopK0sfromAODtracks::SingleTrkCuts(AliAODTrack *trk)
{
  //
  // Single Track Cut to be applied before object creation
  //

  if(trk->GetStatus()&AliESDtrack::kITSpureSA) return kFALSE;
  if(!(trk->GetStatus()&AliESDtrack::kITSin)) return kFALSE;
  if(fProdUseAODFilterBit && !trk->TestFilterMask(BIT(4))) return kFALSE;
  if(fabs(trk->Eta())>fProdTrackEtaRange) return kFALSE;
  if(trk->Pt()<fProdTrackPtMin) return kFALSE;

  if(fUsePID)
    {
      if(fPidHF->GetPidResponse()==0x0){
	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
	AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
	AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
	fPidHF->SetPidResponse(pidResp);
      }

      Int_t isProton=fPidHF->MakeRawPid(trk,4); 
      if(isProton<1) return kFALSE;
    }

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliRDHFCutsLctopK0sfromAODtracks::SingleV0Cuts(AliAODv0 *v0, AliAODVertex *primVert)
{
  //
  // Single V0 Cut to be applied before object creation
  //

  Bool_t onFlyV0 = v0->GetOnFlyStatus(); // on-the-flight V0s
  if ( onFlyV0 && !fUseOnTheFlyV0 ) return kFALSE;

  AliAODTrack *cptrack =  (AliAODTrack*)(v0->GetDaughter(0));
  AliAODTrack *cntrack =  (AliAODTrack*)(v0->GetDaughter(1));
  if(!cptrack || !cntrack) return kFALSE;
  if ( cptrack->Charge() == cntrack->Charge() ) return kFALSE;
  if(!(cptrack->GetStatus() & AliESDtrack::kTPCrefit) ||
     !(cntrack->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
  AliAODVertex *maybeKinkPos = (AliAODVertex*)cptrack->GetProdVertex();
  AliAODVertex *maybeKinkNeg = (AliAODVertex*)cntrack->GetProdVertex();
  if (maybeKinkPos->GetType()==AliAODVertex::kKink || maybeKinkNeg->GetType()==AliAODVertex::kKink) 
    return kFALSE;

  if ( ( ( cptrack->GetTPCClusterInfo(2,1) ) < (Float_t)fProdV0DaughterTPCClusterMin ) || 
       ( ( cntrack->GetTPCClusterInfo(2,1) ) < (Float_t)fProdV0DaughterTPCClusterMin) ) return kFALSE;

  Double_t massK0S = v0->MassK0Short();
  Double_t mk0sPDG   = TDatabasePDG::Instance()->GetParticle(310)->Mass();
  if(fabs(massK0S-mk0sPDG)>fProdV0MassTolK0s) return kFALSE;


  if(TMath::Abs(v0->DcaV0Daughters())>fProdV0DcaDaughtersMax) return kFALSE;
  Double_t posVtx[3] = {0.,0.,0.};
  primVert->GetXYZ(posVtx);
  Double_t cospav0 = v0->CosPointingAngle(posVtx); 
  if(cospav0<fProdV0CosPointingAngleToPrimVtxMin) return kFALSE;
  if(v0->Pt()<fProdV0PtMin) return kFALSE;
  if(fabs(cptrack->Eta())>fProdV0DaughterEtaRange) return kFALSE;
  if(fabs(cntrack->Eta())>fProdV0DaughterEtaRange) return kFALSE;
  if(cptrack->Pt()<fProdV0DaughterPtMin) return kFALSE;
  if(cntrack->Pt()<fProdV0DaughterPtMin) return kFALSE;

  return kTRUE;
}

//________________________________________________________________________
Bool_t AliRDHFCutsLctopK0sfromAODtracks::SelectWithRoughCuts(AliAODv0 *v0, AliAODTrack *part)
{
  //
  // Mass and pT Cut to be applied before object creation
  //

  Double_t mprPDG =  TDatabasePDG::Instance()->GetParticle(2212)->Mass();
  Double_t mLcPDG =  TDatabasePDG::Instance()->GetParticle(4122)->Mass();

  Double_t pxpr_init = part->Px();
  Double_t pypr_init = part->Py();
  Double_t pzpr_init = part->Pz();
  Double_t Epr_init = sqrt(pxpr_init*pxpr_init+pypr_init*pypr_init+pzpr_init*pzpr_init+mprPDG*mprPDG);

  Double_t pxv0_init = v0->Px();
  Double_t pyv0_init = v0->Py();
  Double_t pzv0_init = v0->Pz();
  Double_t massv0_init = v0->MassK0Short();
  Double_t Ev0_init = sqrt(pxv0_init*pxv0_init+pyv0_init*pyv0_init+pzv0_init*pzv0_init+massv0_init*massv0_init);

  Double_t pxlc_init = pxpr_init+pxv0_init;
  Double_t pylc_init = pypr_init+pyv0_init;
  Double_t pzlc_init = pzpr_init+pzv0_init;
  Double_t Elc_init = Epr_init+Ev0_init;
  Double_t lcmass_init = sqrt(Elc_init*Elc_init-pxlc_init*pxlc_init-pylc_init*pylc_init-pzlc_init*pzlc_init);

  if(lcmass_init<mLcPDG-fProdRoughMassTol || lcmass_init>mLcPDG+fProdRoughMassTol) return kFALSE;
  if(sqrt(pxlc_init*pxlc_init+pylc_init*pylc_init)<fProdRoughPtMin) return kFALSE;

  return kTRUE;
}

//________________________________________________________________________
void AliRDHFCutsLctopK0sfromAODtracks::SetCutsArray(Int_t nCuts, Int_t nVars, Int_t nPtBins, Float_t ***cutsArray) 
{
  //
  // store the cuts
  //
  if (nVars!=fnVars) {
    AliError(" wrong number of variables\n");
    return;
  }
  if (nCuts!=fnCuts) {
    AliError(" wrong number of cuts\n");
    return;
  }
  if (nPtBins!=fnPtBins) {
    AliError(" wrong number of pT bin\n");
    return;
  }

  if(!fCutsArray)
    {
      fnTotalCutBins = nCuts * nVars * nPtBins;
      fCutsArray = new Float_t[fnTotalCutBins];
    }

  for(Int_t icuts=0;icuts<nCuts;icuts++){
    for(Int_t iv=0;iv<nVars;iv++){
      for(Int_t ip=0;ip<nPtBins;ip++){
	Int_t gid = GetCutArrayID(icuts,iv,ip);
	fCutsArray[gid] = cutsArray[icuts][iv][ip];
      }
    }
  }
  return;
}

//________________________________________________________________________
void AliRDHFCutsLctopK0sfromAODtracks::SetCutsArray(Int_t nTotBins, Float_t *cutsArray) 
{
  //
  // store the cuts
  //

  if (nTotBins!=fnVars*fnPtBins*fnCuts) {
    AliError(" wrong number of variables\n");
    return;
  }

  if(!fCutsArray)
    {
      fnTotalCutBins = fnCuts * fnVars * fnPtBins;
      fCutsArray = new Float_t[fnTotalCutBins];
    }

  for(Int_t i=0;i<nTotBins;i++){
    fCutsArray[i] = cutsArray[i];
  }

  return;
}

//________________________________________________________________________
void AliRDHFCutsLctopK0sfromAODtracks::SetCutsFromArray(Int_t iCuts) 
{
  //
  // Read Cuts from fCutsArray and update fCutsRD
  //
  Float_t** cutsarray;
  cutsarray=new Float_t*[fnVars];
  for(Int_t ic=0;ic<fnVars;ic++){cutsarray[ic]=new Float_t[fnPtBins];}

  for(Int_t i=0;i<fnVars;i++){
    for(Int_t j=0;j<fnPtBins;j++){
      Int_t gid = GetCutArrayID(iCuts,i,j);
      cutsarray[i][j] = fCutsArray[gid];
    }
  }
  SetCuts(fnVars,fnPtBins,cutsarray);

  for(Int_t ic=0;ic<fnVars;ic++){
    delete[] cutsarray[ic];
  }
  delete[] cutsarray;
}

//________________________________________________________________________
Int_t AliRDHFCutsLctopK0sfromAODtracks::GetCutArrayID(Int_t ic, Int_t iv, Int_t ip) 
{
  //
  // Global index of fCutsArray
  //
  return ic*fnVars*fnPtBins+iv*fnPtBins+ip;
}
Commit	Line	Data
b5c5bbb6	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	/* $Id$ */
	17
	18	/////////////////////////////////////////////////////////////
	19	//
	20	// Class for cuts on AOD reconstructed Lc->p+K0s
	21	//
	22	// Modified by Y.S Watanabe - wyosuke@cns.s.u-tokyo.ac.jp
	23	//
	24	/////////////////////////////////////////////////////////////
	25
	26	#include <Riostream.h>
	27
	28	#include <TDatabasePDG.h>
	29	#include <TMath.h>
	30
	31	#include "AliAnalysisManager.h"
	32	#include "AliInputEventHandler.h"
	33	#include "AliPIDResponse.h"
	34	#include "AliRDHFCutsLctopK0sfromAODtracks.h"
	35	#include "AliAODRecoCascadeHF.h"
	36	#include "AliAODTrack.h"
	37	#include "AliESDtrack.h"
	38	#include "AliESDVertex.h"
	39	#include "AliAODVertex.h"
	40	#include "AliAODv0.h"
	41	#include "AliAODcascade.h"
	42	#include "AliESDv0.h"
	43
	44	using std::cout;
	45	using std::endl;
	46
	47	ClassImp(AliRDHFCutsLctopK0sfromAODtracks)
	48
	49	//--------------------------------------------------------------------------
	50	AliRDHFCutsLctopK0sfromAODtracks::AliRDHFCutsLctopK0sfromAODtracks(const char* name) :
	51	AliRDHFCuts(name),
	52	fPIDStrategy(kNSigmaCuts),
	53	fCombinedPIDThreshold(0.),
	54	fUseOnTheFlyV0(kFALSE),
	55	fProdTrackPtMin(0.3),
	56	fProdTrackEtaRange(0.8),
	57	fProdUseAODFilterBit(kTRUE),
	58	fProdV0MassTolK0s(0.01),
	59	fProdV0PtMin(0.5),
	60	fProdV0CosPointingAngleToPrimVtxMin(0.99),
	61	fProdV0DcaDaughtersMax(1.5),
	62	fProdV0DaughterEtaRange(0.8),
	63	fProdV0DaughterPtMin(0.0),
	64	fProdV0DaughterTPCClusterMin(70),
65	fProdRoughMassTol(0.25),
66	fProdRoughPtMin(0.0),
67	fnCuts(0),
68	fnTotalCutBins(1),
69	fCutsArray(0)
70	{
71	//
72	// Default Constructor
73	//
74
75	const Int_t nvars=7;
76	SetNVars(nvars);
77	TString varNames[nvars]={"Lc inv. mass [GeV/c2]", // 0
78	"Lc pT [GeV/c]", //1
79	"Bachelor pT [GeV/c]", //2
80	"Bachelor d0 [cm]", //3
81	"V0 d0 [cm]", //4
82	"K0s mass [GeV/c2]", //5
83	"Decay Length XY [cm]" //6
84	};
85
86	Bool_t isUpperCut[nvars]={kTRUE, // 0
87	kFALSE, //1
88	kFALSE, //2
89	kTRUE, //3
90	kTRUE, //4
91	kTRUE, //5
92	kFALSE //6
93	};
94	SetVarNames(nvars,varNames,isUpperCut);
95	Bool_t forOpt[nvars]={kFALSE, // 0
96	kFALSE, //1
97	kTRUE, //2
98	kTRUE, //3
99	kTRUE, //4
100	kTRUE, //5
101	kTRUE //6
102	};
103	SetVarsForOpt(nvars,forOpt);
104
105	Float_t limits[2]={0,999999999.};
106	SetPtBins(2,limits);
107	}
108	//--------------------------------------------------------------------------
109	AliRDHFCutsLctopK0sfromAODtracks::AliRDHFCutsLctopK0sfromAODtracks(const AliRDHFCutsLctopK0sfromAODtracks &source) :
110	AliRDHFCuts(source),
111	fPIDStrategy(source.fPIDStrategy),
112	fCombinedPIDThreshold(source.fCombinedPIDThreshold),
113	fUseOnTheFlyV0(source.fUseOnTheFlyV0),
114	fProdTrackPtMin(source.fProdTrackPtMin),
115	fProdTrackEtaRange(source.fProdTrackEtaRange),
116	fProdUseAODFilterBit(source.fProdUseAODFilterBit),
117	fProdV0MassTolK0s(source.fProdV0MassTolK0s),
118	fProdV0PtMin(source.fProdV0PtMin),
119	fProdV0CosPointingAngleToPrimVtxMin(source.fProdV0CosPointingAngleToPrimVtxMin),
120	fProdV0DcaDaughtersMax(source.fProdV0DcaDaughtersMax),
121	fProdV0DaughterEtaRange(source.fProdV0DaughterEtaRange),
122	fProdV0DaughterPtMin(source.fProdV0DaughterPtMin),
123	fProdV0DaughterTPCClusterMin(source.fProdV0DaughterTPCClusterMin),
124	fProdRoughMassTol(source.fProdRoughMassTol),
125	fProdRoughPtMin(source.fProdRoughPtMin),
126	fnCuts(source.fnCuts),
127	fnTotalCutBins(source.fnTotalCutBins),
128	fCutsArray(0)
129	{
130	//
131	// Copy constructor
132	//
133	if(source.fCutsArray) SetCutsArray(source.fnTotalCutBins,source.fCutsArray);
134	}
135	//--------------------------------------------------------------------------
136	AliRDHFCutsLctopK0sfromAODtracks &AliRDHFCutsLctopK0sfromAODtracks::operator=(const AliRDHFCutsLctopK0sfromAODtracks &source)
137	{
138	//
139	// assignment operator
140	//
141
142	if (this != &source) {
143	AliRDHFCuts::operator=(source);
144	}
145
146	fPIDStrategy = source.fPIDStrategy;
147	fCombinedPIDThreshold = source.fCombinedPIDThreshold;
148	fUseOnTheFlyV0 = source.fUseOnTheFlyV0;
149	fProdTrackPtMin = source.fProdTrackPtMin;
150	fProdTrackEtaRange = source.fProdTrackEtaRange;
151	fProdUseAODFilterBit = source.fProdUseAODFilterBit;
152	fProdV0MassTolK0s = source.fProdV0MassTolK0s;
153	fProdV0PtMin = source.fProdV0PtMin;
154	fProdV0CosPointingAngleToPrimVtxMin = source.fProdV0CosPointingAngleToPrimVtxMin;
155	fProdV0DcaDaughtersMax=source.fProdV0DcaDaughtersMax;
156	fProdV0DaughterEtaRange=source.fProdV0DaughterEtaRange;
157	fProdV0DaughterPtMin=source.fProdV0DaughterPtMin;
158	fProdV0DaughterTPCClusterMin=source.fProdV0DaughterTPCClusterMin;
159	fProdRoughMassTol = source.fProdRoughMassTol;
160	fProdRoughPtMin = source.fProdRoughPtMin;
161
162	if(source.fnCuts) fnCuts = source.fnCuts;
163	if(source.fnTotalCutBins) fnTotalCutBins = source.fnCuts;
164	if(source.fCutsArray) SetCutsArray(source.fnTotalCutBins,source.fCutsArray);
165
166	return *this;
167	}
168
169	//---------------------------------------------------------------------------
170	AliRDHFCutsLctopK0sfromAODtracks::~AliRDHFCutsLctopK0sfromAODtracks() {
171	//
172	// Default Destructor
173	//
174
175	if(fCutsArray) delete[] fCutsArray;
176	}
177
178	//---------------------------------------------------------------------------
179	void AliRDHFCutsLctopK0sfromAODtracks::GetCutVarsForOpt(AliAODRecoDecayHF d,Float_t vars,Int_t nvars,Int_t *pdgdaughters) {
180	//
181	// Fills in vars the values of the variables
182	//
183
184	if (pdgdaughters[0]==-9999) return; // dummy
185
186	AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)d;
187	if(!dd){
188	AliError("No AliAODRecoCascadeHF object found\n");
189	return;
190	}
191
192	if (nvars!=fnVarsForOpt) {
193	AliError("Cut object seems to have the wrong number of variables\n");
194	return;
195	}
196
197	//Double_t ptD=d->Pt();
198	//Int_t ptbin=PtBin(ptD);
199	Int_t iter=-1;
200
201	if(fVarsForOpt[0]){
202	iter++;
203	Double_t mlcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
204	vars[iter]= TMath::Abs(dd->InvMassLctoK0sP()-mlcPDG) ;
205	}
206	if(fVarsForOpt[1]){
207	iter++;
208	vars[iter]= dd->Pt();
209	}
210	if(fVarsForOpt[2]){
211	iter++;
212	AliAODTrack *part = dd->GetBachelor();
213	vars[iter]= part->Pt();
214	}
215	if(fVarsForOpt[3]){
216	iter++;
217	vars[iter]= dd->Getd0Prong(0);
218	}
219	if(fVarsForOpt[4]){
220	iter++;
221	vars[iter]= dd->Getd0Prong(1);
222	}
223	if(fVarsForOpt[5]){
224	iter++;
225	AliAODv0 *v0 = dd->Getv0();
226	vars[iter]= v0->MassK0Short();
227	}
228	if(fVarsForOpt[6]){
229	iter++;
230	vars[iter]= dd->DecayLengthXY();
231	}
232
233	return;
234	}
235
236	//---------------------------------------------------------------------------
237	Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelected(TObject* obj, Int_t selectionLevel) {
238	//
239	// Apply selection
240	//
241
242	if (!fCutsRD) {
243	AliFatal("Cut matrix not inizialized. Exit...");
244	return 0;
245	}
246
247	AliAODRecoCascadeHF* d=(AliAODRecoCascadeHF*)obj;
248	if(!d){
249	AliDebug(2,"AliAODRecoCascadeHF null");
250	return 0;
251	}
252
253	Double_t ptD=d->Pt();
254	if(ptD<fMinPtCand) return 0;
255	if(ptD>fMaxPtCand) return 0;
256
257	if (selectionLevel==AliRDHFCuts::kAll \|\|
258	selectionLevel==AliRDHFCuts::kTracks) {
259	//Performed in production stage
260	}
261
262	Int_t returnvalueCuts=1;
263	// selection on candidate
264	if (selectionLevel==AliRDHFCuts::kAll \|\|
265	selectionLevel==AliRDHFCuts::kCandidate) {
266
267	Double_t pt=d->Pt();
268	Int_t ptbin=PtBin(pt);
269	if (ptbin==-1) {
270	return 0;
271	}
272	Bool_t okcand=kTRUE;
273
274	Double_t mlcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
275	Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
276	AliAODTrack *part = d->GetBachelor();
277	AliAODv0 *v0 = d->Getv0();
278
279	if(TMath::Abs(d->InvMassLctoK0sP()-mlcPDG) > fCutsRD[GetGlobalIndex(0,ptbin)])
280	{
281	okcand = kFALSE;
282	}
283	if(d->Pt() < fCutsRD[GetGlobalIndex(1,ptbin)])
284	{
285	okcand = kFALSE;
286	}
287	if(part->Pt() < fCutsRD[GetGlobalIndex(2,ptbin)])
288	{
289	okcand = kFALSE;
290	}
291	if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(3,ptbin)])
292	{
293	okcand = kFALSE;
294	}
295	if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(4,ptbin)])
296	{
297	okcand = kFALSE;
298	}
299	if(TMath::Abs(v0->MassK0Short()-mk0sPDG) > fCutsRD[GetGlobalIndex(5,ptbin)])
300	{
301	okcand = kFALSE;
302	}
303	if(d->DecayLengthXY() < fCutsRD[GetGlobalIndex(6,ptbin)])
304	{
305	okcand = kFALSE;
306	}
307
308	if(!okcand) return 0;
309	returnvalueCuts = 1;
310	}
311
312	Int_t returnvaluePID=1;
313	if(selectionLevel==AliRDHFCuts::kAll \|\|
314	selectionLevel==AliRDHFCuts::kCandidate\|\|
315	selectionLevel==AliRDHFCuts::kPID) {
316
317	switch(fPIDStrategy){
318	case kNSigmaCuts:
319	returnvaluePID = IsSelectedPID(d);
320	break;
321	case kCombinedCuts:
322	returnvaluePID = IsSelectedCombinedPID(d);
323	break;
324	}
325	}
326
327	Int_t returnvalue = 0;
328	if(returnvalueCuts==1 && returnvaluePID==1) returnvalue=1;
329
330	return returnvalue;
331	}
332
333	//---------------------------------------------------------------------------
334	Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelectedPID(AliAODRecoDecayHF* obj)
335	{
336	//
337	// IsSelectedPID
338	//
339
340	if(!fUsePID \|\| !obj) return 1;
341
342	AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)obj;
343	AliAODTrack *part = dd->GetBachelor();
344
345	Int_t returnvalue=1;
346
347	if(fPidHF->GetPidResponse()==0x0){
348	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
349	AliInputEventHandler inputHandler=(AliInputEventHandler)mgr->GetInputEventHandler();
350	AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
351	fPidHF->SetPidResponse(pidResp);
352	}
353
354	Int_t isProton=fPidHF->MakeRawPid(part,4);
355	if(isProton<1) returnvalue = 0;
356
357	return returnvalue;
358	}
359
360	//---------------------------------------------------------------------------
361	Int_t AliRDHFCutsLctopK0sfromAODtracks::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
362	//
363	// IsSelectedCombinedPID
364	//
365
366	if(!fUsePID \|\| !obj) {return 1;}
367
368	AliAODRecoCascadeHF* dd=(AliAODRecoCascadeHF*)obj;
369	AliAODTrack *part = dd->GetBachelor();
370	if(!part) return 0;
371
372	Int_t returnvalue=1;
373	Double_t probProton = GetProtonProbabilityTPCTOF(part);
374	if(probProton<fCombinedPIDThreshold) returnvalue = 0;
375
376	return returnvalue;
377	}
378
379	//________________________________________________________________________
380	Double_t AliRDHFCutsLctopK0sfromAODtracks::GetProtonProbabilityTPCTOF(AliAODTrack* trk)
381	{
382	//
383	// Get Proton probability
384	//
385	fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC+AliPIDResponse::kDetTOF);
386	Double_t prob1[AliPID::kSPECIES];
387	UInt_t detUsed1 = fPidHF->GetPidCombined()->ComputeProbabilities(trk, fPidHF->GetPidResponse(), prob1);
388	if (detUsed1 != (UInt_t)fPidHF->GetPidCombined()->GetDetectorMask() )
389	{
390	fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
391	detUsed1 = fPidHF->GetPidCombined()->ComputeProbabilities(trk, fPidHF->GetPidResponse(), prob1);
392	}
393	return prob1[AliPID::kProton];
394	}
395
396	//________________________________________________________________________
397	Bool_t AliRDHFCutsLctopK0sfromAODtracks::SingleTrkCuts(AliAODTrack *trk)
398	{
399	//
400	// Single Track Cut to be applied before object creation
401	//
402
403	if(trk->GetStatus()&AliESDtrack::kITSpureSA) return kFALSE;
404	if(!(trk->GetStatus()&AliESDtrack::kITSin)) return kFALSE;
405	if(fProdUseAODFilterBit && !trk->TestFilterMask(BIT(4))) return kFALSE;
406	if(fabs(trk->Eta())>fProdTrackEtaRange) return kFALSE;
407	if(trk->Pt()<fProdTrackPtMin) return kFALSE;
408
409	if(fUsePID)
410	{
411	if(fPidHF->GetPidResponse()==0x0){
412	AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
413	AliInputEventHandler inputHandler=(AliInputEventHandler)mgr->GetInputEventHandler();
414	AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
415	fPidHF->SetPidResponse(pidResp);
416	}
417
418	Int_t isProton=fPidHF->MakeRawPid(trk,4);
419	if(isProton<1) return kFALSE;
420	}
421
422	return kTRUE;
423	}
424
425	//________________________________________________________________________
426	Bool_t AliRDHFCutsLctopK0sfromAODtracks::SingleV0Cuts(AliAODv0 v0, AliAODVertex primVert)
427	{
428	//
429	// Single V0 Cut to be applied before object creation
430	//
431
432	Bool_t onFlyV0 = v0->GetOnFlyStatus(); // on-the-flight V0s
433	if ( onFlyV0 && !fUseOnTheFlyV0 ) return kFALSE;
434
435	AliAODTrack cptrack = (AliAODTrack)(v0->GetDaughter(0));
436	AliAODTrack cntrack = (AliAODTrack)(v0->GetDaughter(1));
437	if(!cptrack \|\| !cntrack) return kFALSE;
438	if ( cptrack->Charge() == cntrack->Charge() ) return kFALSE;
439	if(!(cptrack->GetStatus() & AliESDtrack::kTPCrefit) \|\|
440	!(cntrack->GetStatus() & AliESDtrack::kTPCrefit)) return kFALSE;
441	AliAODVertex maybeKinkPos = (AliAODVertex)cptrack->GetProdVertex();
442	AliAODVertex maybeKinkNeg = (AliAODVertex)cntrack->GetProdVertex();
443	if (maybeKinkPos->GetType()==AliAODVertex::kKink \|\| maybeKinkNeg->GetType()==AliAODVertex::kKink)
444	return kFALSE;
445
446	if ( ( ( cptrack->GetTPCClusterInfo(2,1) ) < (Float_t)fProdV0DaughterTPCClusterMin ) \|\|
447	( ( cntrack->GetTPCClusterInfo(2,1) ) < (Float_t)fProdV0DaughterTPCClusterMin) ) return kFALSE;
448
449	Double_t massK0S = v0->MassK0Short();
450	Double_t mk0sPDG = TDatabasePDG::Instance()->GetParticle(310)->Mass();
451	if(fabs(massK0S-mk0sPDG)>fProdV0MassTolK0s) return kFALSE;
452
453
454	if(TMath::Abs(v0->DcaV0Daughters())>fProdV0DcaDaughtersMax) return kFALSE;
455	Double_t posVtx[3] = {0.,0.,0.};
456	primVert->GetXYZ(posVtx);
457	Double_t cospav0 = v0->CosPointingAngle(posVtx);
458	if(cospav0<fProdV0CosPointingAngleToPrimVtxMin) return kFALSE;
459	if(v0->Pt()<fProdV0PtMin) return kFALSE;
460	if(fabs(cptrack->Eta())>fProdV0DaughterEtaRange) return kFALSE;
461	if(fabs(cntrack->Eta())>fProdV0DaughterEtaRange) return kFALSE;
462	if(cptrack->Pt()<fProdV0DaughterPtMin) return kFALSE;
463	if(cntrack->Pt()<fProdV0DaughterPtMin) return kFALSE;
464
465	return kTRUE;
466	}
467
468	//________________________________________________________________________
469	Bool_t AliRDHFCutsLctopK0sfromAODtracks::SelectWithRoughCuts(AliAODv0 v0, AliAODTrack part)
470	{
471	//
472	// Mass and pT Cut to be applied before object creation
473	//
474
475	Double_t mprPDG = TDatabasePDG::Instance()->GetParticle(2212)->Mass();
476	Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
477
478	Double_t pxpr_init = part->Px();
479	Double_t pypr_init = part->Py();
480	Double_t pzpr_init = part->Pz();
481	Double_t Epr_init = sqrt(pxpr_initpxpr_init+pypr_initpypr_init+pzpr_initpzpr_init+mprPDGmprPDG);
482
483	Double_t pxv0_init = v0->Px();
484	Double_t pyv0_init = v0->Py();
485	Double_t pzv0_init = v0->Pz();
486	Double_t massv0_init = v0->MassK0Short();
487	Double_t Ev0_init = sqrt(pxv0_initpxv0_init+pyv0_initpyv0_init+pzv0_initpzv0_init+massv0_initmassv0_init);
488
489	Double_t pxlc_init = pxpr_init+pxv0_init;
490	Double_t pylc_init = pypr_init+pyv0_init;
491	Double_t pzlc_init = pzpr_init+pzv0_init;
492	Double_t Elc_init = Epr_init+Ev0_init;
493	Double_t lcmass_init = sqrt(Elc_initElc_init-pxlc_initpxlc_init-pylc_initpylc_init-pzlc_initpzlc_init);
494
495	if(lcmass_init<mLcPDG-fProdRoughMassTol \|\| lcmass_init>mLcPDG+fProdRoughMassTol) return kFALSE;
496	if(sqrt(pxlc_initpxlc_init+pylc_initpylc_init)<fProdRoughPtMin) return kFALSE;
497
498	return kTRUE;
499	}
500
501	//________________________________________________________________________
502	void AliRDHFCutsLctopK0sfromAODtracks::SetCutsArray(Int_t nCuts, Int_t nVars, Int_t nPtBins, Float_t ***cutsArray)
503	{
504	//
505	// store the cuts
506	//
507	if (nVars!=fnVars) {
508	AliError(" wrong number of variables\n");
509	return;
510	}
511	if (nCuts!=fnCuts) {
512	AliError(" wrong number of cuts\n");
513	return;
514	}
515	if (nPtBins!=fnPtBins) {
516	AliError(" wrong number of pT bin\n");
517	return;
518	}
519
520	if(!fCutsArray)
521	{
522	fnTotalCutBins = nCuts * nVars * nPtBins;
523	fCutsArray = new Float_t[fnTotalCutBins];
524	}
525
526	for(Int_t icuts=0;icuts<nCuts;icuts++){
527	for(Int_t iv=0;iv<nVars;iv++){
528	for(Int_t ip=0;ip<nPtBins;ip++){
529	Int_t gid = GetCutArrayID(icuts,iv,ip);
530	fCutsArray[gid] = cutsArray[icuts][iv][ip];
531	}
532	}
533	}
534	return;
535	}
536
537	//________________________________________________________________________
538	void AliRDHFCutsLctopK0sfromAODtracks::SetCutsArray(Int_t nTotBins, Float_t *cutsArray)
539	{
540	//
541	// store the cuts
542	//
543
544	if (nTotBins!=fnVarsfnPtBinsfnCuts) {
545	AliError(" wrong number of variables\n");
546	return;
547	}
548
549	if(!fCutsArray)
550	{
551	fnTotalCutBins = fnCuts * fnVars * fnPtBins;
552	fCutsArray = new Float_t[fnTotalCutBins];
553	}
554
555	for(Int_t i=0;i<nTotBins;i++){
556	fCutsArray[i] = cutsArray[i];
557	}
558
559	return;
560	}
561
562	//________________________________________________________________________
563	void AliRDHFCutsLctopK0sfromAODtracks::SetCutsFromArray(Int_t iCuts)
564	{
565	//
566	// Read Cuts from fCutsArray and update fCutsRD
567	//
568	Float_t** cutsarray;
569	cutsarray=new Float_t*[fnVars];
570	for(Int_t ic=0;ic<fnVars;ic++){cutsarray[ic]=new Float_t[fnPtBins];}
571
572	for(Int_t i=0;i<fnVars;i++){
573	for(Int_t j=0;j<fnPtBins;j++){
574	Int_t gid = GetCutArrayID(iCuts,i,j);
575	cutsarray[i][j] = fCutsArray[gid];
576	}
577	}
578	SetCuts(fnVars,fnPtBins,cutsarray);
579
580	for(Int_t ic=0;ic<fnVars;ic++){
581	delete[] cutsarray[ic];
582	}
583	delete[] cutsarray;
584	}
585
586	//________________________________________________________________________
587	Int_t AliRDHFCutsLctopK0sfromAODtracks::GetCutArrayID(Int_t ic, Int_t iv, Int_t ip)
588	{
589	//
590	// Global index of fCutsArray
591	//
592	return icfnVarsfnPtBins+iv*fnPtBins+ip;
593	}