-/**************************************************************************\r
- * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-\r
-/* $Id$ */\r
-\r
-/////////////////////////////////////////////////////////////\r
-//\r
-// Class for cuts on AOD reconstructed D0->Kpi\r
-//\r
-// Author: A.Dainese, andrea.dainese@pd.infn.it\r
-/////////////////////////////////////////////////////////////\r
-\r
-#include <TDatabasePDG.h>\r
-#include <Riostream.h>\r
-\r
-#include "AliRDHFCutsD0toKpi.h"\r
-#include "AliAODRecoDecayHF2Prong.h"\r
-#include "AliAODTrack.h"\r
-#include "AliESDtrack.h"\r
-#include "AliAODPid.h"\r
-#include "AliTPCPIDResponse.h"\r
-#include "AliAODVertex.h"\r
-#include "AliKFVertex.h"\r
-#include "AliKFParticle.h"\r
-\r
-ClassImp(AliRDHFCutsD0toKpi)\r
-\r
-//--------------------------------------------------------------------------\r
-AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) : \r
-AliRDHFCuts(name),\r
-fUseSpecialCuts(kFALSE),\r
-fLowPt(kTRUE),\r
-fDefaultPID(kFALSE),\r
-fUseKF(kFALSE),\r
-fPtLowPID(2.),\r
-fPtMaxSpecialCuts(9999.)\r
-{\r
- //\r
- // Default Constructor\r
- //\r
- Int_t nvars=11;\r
- SetNVars(nvars);\r
- TString varNames[11]={"inv. mass [GeV]", \r
- "dca [cm]",\r
- "cosThetaStar", \r
- "pTK [GeV/c]",\r
- "pTPi [GeV/c]",\r
- "d0K [cm]",\r
- "d0Pi [cm]",\r
- "d0d0 [cm^2]",\r
- "cosThetaPoint",\r
- "|cosThetaPointXY|", \r
- "NormDecayLenghtXY"};\r
- Bool_t isUpperCut[11]={kTRUE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE, \r
- kFALSE};\r
- SetVarNames(nvars,varNames,isUpperCut);\r
- Bool_t forOpt[11]={kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE,\r
- kFALSE,\r
- kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE};\r
- SetVarsForOpt(4,forOpt);\r
- Float_t limits[2]={0,999999999.};\r
- SetPtBins(2,limits);\r
-\r
-}\r
-//--------------------------------------------------------------------------\r
-AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :\r
- AliRDHFCuts(source), \r
- fUseSpecialCuts(source.fUseSpecialCuts),\r
- fLowPt(source.fLowPt),\r
- fDefaultPID(source.fDefaultPID),\r
- fUseKF(source.fUseKF),\r
- fPtLowPID(source.fPtLowPID),\r
- fPtMaxSpecialCuts(source.fPtMaxSpecialCuts)\r
-{\r
- //\r
- // Copy constructor\r
- //\r
-\r
-}\r
-//--------------------------------------------------------------------------\r
-AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)\r
-{\r
- //\r
- // assignment operator\r
- //\r
- if(&source == this) return *this;\r
-\r
- AliRDHFCuts::operator=(source); \r
- fUseSpecialCuts=source.fUseSpecialCuts;\r
- fLowPt=source.fLowPt;\r
- fDefaultPID=source.fDefaultPID;\r
- fUseKF=source.fUseKF;\r
- fPtLowPID=source.fPtLowPID;\r
- fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;\r
-\r
- return *this;\r
-}\r
-\r
-\r
-//---------------------------------------------------------------------------\r
-void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {\r
- // \r
- // Fills in vars the values of the variables \r
- //\r
-\r
- if(nvars!=fnVarsForOpt) {\r
- printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");\r
- return;\r
- }\r
-\r
- AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;\r
- \r
- //recalculate vertex w/o daughters\r
- Bool_t cleanvtx=kFALSE;\r
- AliAODVertex *origownvtx=0x0;\r
- if(fRemoveDaughtersFromPrimary) {\r
- if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());\r
- cleanvtx=kTRUE;\r
- if(!RecalcOwnPrimaryVtx(dd,aod)) {\r
- CleanOwnPrimaryVtx(dd,aod,origownvtx);\r
- cleanvtx=kFALSE;\r
- }\r
- }\r
-\r
- Int_t iter=-1;\r
- if(fVarsForOpt[0]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==211) {\r
- vars[iter]=dd->InvMassD0();\r
- } else {\r
- vars[iter]=dd->InvMassD0bar();\r
- }\r
- }\r
- if(fVarsForOpt[1]){\r
- iter++;\r
- vars[iter]=dd->GetDCA();\r
- }\r
- if(fVarsForOpt[2]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==211) {\r
- vars[iter] = dd->CosThetaStarD0();\r
- } else {\r
- vars[iter] = dd->CosThetaStarD0bar();\r
- }\r
- }\r
- if(fVarsForOpt[3]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==321) {\r
- vars[iter]=dd->PtProng(0);\r
- }\r
- else{\r
- vars[iter]=dd->PtProng(1);\r
- }\r
- }\r
- if(fVarsForOpt[4]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==211) {\r
- vars[iter]=dd->PtProng(0);\r
- }\r
- else{\r
- vars[iter]=dd->PtProng(1);\r
- }\r
- }\r
- if(fVarsForOpt[5]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==321) {\r
- vars[iter]=dd->Getd0Prong(0);\r
- }\r
- else{\r
- vars[iter]=dd->Getd0Prong(1);\r
- }\r
- }\r
- if(fVarsForOpt[6]){\r
- iter++;\r
- if(TMath::Abs(pdgdaughters[0])==211) {\r
- vars[iter]=dd->Getd0Prong(0);\r
- }\r
- else{\r
- vars[iter]=dd->Getd0Prong(1);\r
- }\r
- }\r
- if(fVarsForOpt[7]){\r
- iter++;\r
- vars[iter]= dd->Prodd0d0();\r
- }\r
- if(fVarsForOpt[8]){\r
- iter++;\r
- vars[iter]=dd->CosPointingAngle();\r
- }\r
- \r
- if(fVarsForOpt[9]){\r
- iter++;\r
- vars[iter]=TMath::Abs(dd->CosPointingAngleXY());\r
- }\r
- \r
- if(fVarsForOpt[10]){\r
- iter++;\r
- vars[iter]=dd->NormalizedDecayLengthXY();\r
- }\r
-\r
- if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);\r
-\r
- return;\r
-}\r
-//---------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {\r
- //\r
- // Apply selection\r
- //\r
- \r
-\r
- fIsSelectedCuts=0;\r
- fIsSelectedPID=0;\r
-\r
- if(!fCutsRD){\r
- cout<<"Cut matrice not inizialized. Exit..."<<endl;\r
- return 0;\r
- }\r
- //PrintAll();\r
- AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;\r
-\r
- if(!d){\r
- cout<<"AliAODRecoDecayHF2Prong null"<<endl;\r
- return 0;\r
- }\r
-\r
- if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;\r
-\r
- Double_t ptD=d->Pt();\r
- if(ptD<fMinPtCand) return 0;\r
- if(ptD>fMaxPtCand) return 0;\r
-\r
- if(d->HasBadDaughters()) return 0;\r
-\r
- // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both\r
- Int_t returnvaluePID=3;\r
- Int_t returnvalueCuts=3;\r
-\r
- // selection on candidate\r
- if(selectionLevel==AliRDHFCuts::kAll || \r
- selectionLevel==AliRDHFCuts::kCandidate) {\r
-\r
- if(!fUseKF) {\r
-\r
- //recalculate vertex w/o daughters\r
- AliAODVertex *origownvtx=0x0;\r
- if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {\r
- if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());\r
- if(!RecalcOwnPrimaryVtx(d,aod)) { \r
- CleanOwnPrimaryVtx(d,aod,origownvtx);\r
- return 0;\r
- }\r
- }\r
-\r
- if(fUseMCVertex) {\r
- if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());\r
- if(!SetMCPrimaryVtx(d,aod)) {\r
- CleanOwnPrimaryVtx(d,aod,origownvtx);\r
- return 0;\r
- }\r
- }\r
- \r
- Double_t pt=d->Pt();\r
- \r
- Int_t okD0=0,okD0bar=0;\r
- \r
- Int_t ptbin=PtBin(pt);\r
- if (ptbin==-1) {\r
- CleanOwnPrimaryVtx(d,aod,origownvtx);\r
- return 0;\r
- }\r
-\r
- Double_t mD0,mD0bar,ctsD0,ctsD0bar;\r
- okD0=1; okD0bar=1;\r
- \r
- Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
-\r
- d->InvMassD0(mD0,mD0bar);\r
- if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;\r
- if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
-\r
- if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- \r
- if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;\r
- if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- \r
- if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || \r
- TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;\r
- if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||\r
- TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- \r
- d->CosThetaStarD0(ctsD0,ctsD0bar);\r
- if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; \r
- if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
-\r
- \r
- if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- Double_t normalDecayLengXY=d->NormalizedDecayLengthXY();\r
- if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
- \r
- if (returnvalueCuts!=0) {\r
- if (okD0) returnvalueCuts=1; //cuts passed as D0\r
- if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar\r
- if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar\r
- }\r
-\r
- // call special cuts\r
- Int_t special=1;\r
- if(fUseSpecialCuts && (pt<fPtMaxSpecialCuts)) special=IsSelectedSpecialCuts(d);\r
- if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}\r
-\r
- // unset recalculated primary vertex when not needed any more\r
- CleanOwnPrimaryVtx(d,aod,origownvtx);\r
-\r
- } else {\r
- // go to selection with Kalman vertexing, if requested\r
- returnvalueCuts = IsSelectedKF(d,aod);\r
- }\r
-\r
- fIsSelectedCuts=returnvalueCuts;\r
- if(!returnvalueCuts) return 0;\r
- }\r
-\r
-\r
- // selection on PID \r
- if(selectionLevel==AliRDHFCuts::kAll || \r
- selectionLevel==AliRDHFCuts::kCandidate ||\r
- selectionLevel==AliRDHFCuts::kPID) {\r
- returnvaluePID = IsSelectedPID(d);\r
- fIsSelectedPID=returnvaluePID;\r
- if(!returnvaluePID) return 0;\r
- }\r
-\r
- Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);\r
-\r
- if(!returnvalueComb) return 0;\r
-\r
- // selection on daughter tracks \r
- if(selectionLevel==AliRDHFCuts::kAll || \r
- selectionLevel==AliRDHFCuts::kTracks) {\r
- if(!AreDaughtersSelected(d)) return 0;\r
- }\r
- \r
- // cout<<"Pid = "<<returnvaluePID<<endl;\r
- return returnvalueComb;\r
-}\r
-\r
-//------------------------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,\r
- AliAODEvent* aod) const {\r
- //\r
- // Apply selection using KF-vertexing\r
- //\r
- \r
- AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);\r
- AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1); \r
- \r
- if(!track0 || !track1) {\r
- cout<<"one or two D0 daughters missing!"<<endl;\r
- return 0;\r
- }\r
-\r
- // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both\r
- Int_t returnvalueCuts=3;\r
- \r
- // candidate selection with AliKF\r
- AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field\r
- \r
- Int_t okD0=0,okD0bar=0;\r
- okD0=1; okD0bar=1;\r
- \r
- // convert tracks into AliKFParticles\r
- \r
- AliKFParticle negPiKF(*track1,-211); // neg pion kandidate\r
- AliKFParticle negKKF(*track1,-321); // neg kaon kandidate\r
- AliKFParticle posPiKF(*track0,211); // pos pion kandidate\r
- AliKFParticle posKKF(*track0,321); // pos kaon kandidate\r
- \r
- // build D0 candidates\r
- \r
- AliKFParticle d0c(negKKF,posPiKF); // D0 candidate\r
- AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate\r
- \r
- // create kf primary vertices\r
- \r
- AliAODVertex *vtx1 = aod->GetPrimaryVertex();\r
- AliKFVertex primVtx(*vtx1); \r
- AliKFVertex aprimVtx(*vtx1);\r
- \r
- if(primVtx.GetNContributors()<=0) okD0 = 0;\r
- if(aprimVtx.GetNContributors()<=0) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- // calculate mass\r
- \r
- Double_t d0mass = d0c.GetMass();\r
- Double_t ad0mass = ad0c.GetMass();\r
- \r
- // calculate P of D0 and D0bar\r
- Double_t d0P = d0c.GetP();\r
- Double_t d0Px = d0c.GetPx();\r
- Double_t d0Py = d0c.GetPy();\r
- Double_t d0Pz = d0c.GetPz();\r
- Double_t ad0P = ad0c.GetP(); \r
- Double_t ad0Px = ad0c.GetPx();\r
- Double_t ad0Py = ad0c.GetPy();\r
- Double_t ad0Pz = ad0c.GetPz();\r
- \r
- //calculate Pt of D0 and D0bar\r
- \r
- Double_t pt=d0c.GetPt(); \r
- Double_t apt=ad0c.GetPt();\r
- \r
- // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates\r
- \r
- if(track0->GetUsedForPrimVtxFit()) {\r
- primVtx -= posPiKF; \r
- aprimVtx -= posKKF;\r
- }\r
- \r
- if(track1->GetUsedForPrimVtxFit()) { \r
- primVtx -= negKKF; \r
- aprimVtx -= negPiKF;\r
- }\r
- \r
- primVtx += d0c;\r
- aprimVtx += ad0c;\r
- \r
- if(primVtx.GetNContributors()<=0) okD0 = 0;\r
- if(aprimVtx.GetNContributors()<=0) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- //calculate cut variables\r
- \r
- // calculate impact params of daughters w.r.t recalculated vertices\r
- \r
- Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);\r
- Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);\r
- Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);\r
- Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);\r
- \r
- // calculate Product of Impact Params\r
- \r
- Double_t prodParam = impactPi*impactKa;\r
- Double_t aprodParam = aimpactPi*aimpactKa;\r
- \r
- // calculate cosine of pointing angles\r
- \r
- TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());\r
- TVector3 fline(d0c.GetX()-primVtx.GetX(),\r
- d0c.GetY()-primVtx.GetY(),\r
- d0c.GetZ()-primVtx.GetZ());\r
- Double_t pta = mom.Angle(fline);\r
- Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate\r
- \r
- TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());\r
- TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),\r
- ad0c.GetY()-aprimVtx.GetY(),\r
- ad0c.GetZ()-aprimVtx.GetZ());\r
- Double_t apta = amom.Angle(afline);\r
- Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate\r
- \r
- // calculate P of Pions at Decay Position of D0 and D0bar candidates\r
- negKKF.TransportToParticle(d0c);\r
- posPiKF.TransportToParticle(d0c);\r
- posKKF.TransportToParticle(ad0c);\r
- negPiKF.TransportToParticle(ad0c);\r
- \r
- Double_t pxPi = posPiKF.GetPx();\r
- Double_t pyPi = posPiKF.GetPy();\r
- Double_t pzPi = posPiKF.GetPz();\r
- Double_t ptPi = posPiKF.GetPt();\r
- \r
- Double_t apxPi = negPiKF.GetPx();\r
- Double_t apyPi = negPiKF.GetPy();\r
- Double_t apzPi = negPiKF.GetPz();\r
- Double_t aptPi = negPiKF.GetPt();\r
- \r
- // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates\r
- \r
- Double_t ptK = negKKF.GetPt();\r
- Double_t aptK = posKKF.GetPt();\r
- \r
- //calculate cos(thetastar)\r
- Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();\r
- Double_t massp[2];\r
- massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();\r
- massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();\r
- Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)\r
- -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);\r
- \r
- // cos(thetastar) for D0 and Pion\r
- \r
- Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);\r
- Double_t beta = d0P/d0E;\r
- Double_t gamma = d0E/massvtx;\r
- TVector3 momPi(pxPi,pyPi,pzPi);\r
- TVector3 momTot(d0Px,d0Py,d0Pz);\r
- Double_t q1 = momPi.Dot(momTot)/momTot.Mag();\r
- Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar; \r
- \r
- // cos(thetastar) for D0bar and Pion \r
- \r
- Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);\r
- Double_t abeta = ad0P/ad0E;\r
- Double_t agamma = ad0E/massvtx;\r
- TVector3 amomPi(apxPi,apyPi,apzPi);\r
- TVector3 amomTot(ad0Px,ad0Py,ad0Pz);\r
- Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();\r
- Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar; \r
- \r
- // calculate reduced Chi2 for the full D0 fit\r
- d0c.SetProductionVertex(primVtx);\r
- ad0c.SetProductionVertex(aprimVtx);\r
- negKKF.SetProductionVertex(d0c);\r
- posPiKF.SetProductionVertex(d0c);\r
- posKKF.SetProductionVertex(ad0c);\r
- negPiKF.SetProductionVertex(ad0c);\r
- d0c.TransportToProductionVertex();\r
- ad0c.TransportToProductionVertex();\r
- \r
- // calculate the decay length\r
- Double_t decayLengthD0 = d0c.GetDecayLength();\r
- Double_t adecayLengthD0 = ad0c.GetDecayLength();\r
- \r
- Double_t chi2D0 = 50.;\r
- if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {\r
- chi2D0 = d0c.GetChi2()/d0c.GetNDF();\r
- }\r
- \r
- Double_t achi2D0 = 50.;\r
- if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {\r
- achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();\r
- }\r
- \r
- // Get the Pt-bins\r
- Int_t ptbin=PtBin(pt);\r
- Int_t aptbin=PtBin(apt);\r
-\r
- if(ptbin < 0) okD0 = 0;\r
- if(aptbin < 0) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;\r
- if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- \r
- if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] || \r
- TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;\r
- \r
- if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||\r
- TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;\r
- \r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- // for the moment via the standard method due to bug in AliKF\r
- if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;\r
- if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- \r
- if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;\r
- if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- \r
- if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; \r
- if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;\r
- if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0; \r
- if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0; \r
- if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0; \r
- if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;\r
- if(!okD0 && !okD0bar) returnvalueCuts=0;\r
- \r
- if(returnvalueCuts!=0) {\r
- if(okD0) returnvalueCuts=1; //cuts passed as D0\r
- if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar\r
- if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar\r
- }\r
-\r
- return returnvalueCuts; \r
-}\r
-\r
-//---------------------------------------------------------------------------\r
-\r
-Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const\r
-{\r
- //\r
- // Checking if D0 is in fiducial acceptance region \r
- //\r
-\r
- if(pt > 5.) {\r
- // applying cut for pt > 5 GeV\r
- AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt)); \r
- if (TMath::Abs(y) > 0.8){\r
- return kFALSE;\r
- }\r
- } else {\r
- // appliying smooth cut for pt < 5 GeV\r
- Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5; \r
- Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5; \r
- AliDebug(2,Form("pt of D0 = %f (< 5), cutting according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY)); \r
- if (y < minFiducialY || y > maxFiducialY){\r
- return kFALSE;\r
- }\r
- }\r
-\r
- return kTRUE;\r
-}\r
-//---------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d) \r
-{\r
- // ############################################################\r
- //\r
- // Apply PID selection\r
- //\r
- //\r
- // ############################################################\r
-\r
- if(!fUsePID) return 3;\r
- if(fDefaultPID) return IsSelectedPIDdefault(d);\r
- fWhyRejection=0;\r
- Int_t isD0D0barPID[2]={1,2};\r
- Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value; \r
- // same for prong 2\r
- // values convention -1 = discarded \r
- // 0 = not identified (but compatible) || No PID (->hasPID flag)\r
- // 1 = identified\r
- // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both \r
- // Initial hypothesis: unknwon (but compatible) \r
- combinedPID[0][0]=0; // prima figlia, Kaon\r
- combinedPID[0][1]=0; // prima figlia, pione\r
- combinedPID[1][0]=0; // seconda figlia, Kaon\r
- combinedPID[1][1]=0; // seconda figlia, pion\r
- \r
- Bool_t checkPIDInfo[2]={kTRUE,kTRUE};\r
- Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};\r
- for(Int_t daught=0;daught<2;daught++){\r
- //Loop con prongs\r
- AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);\r
- if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0; \r
- \r
- if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {\r
- checkPIDInfo[daught]=kFALSE; \r
- continue;\r
- }\r
-\r
- // identify kaon\r
- combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);\r
-\r
- // identify pion\r
-\r
- if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {\r
- combinedPID[daught][1]=0;\r
- }else{\r
- fPidHF->SetTOF(kFALSE);\r
- combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);\r
- fPidHF->SetTOF(kTRUE);\r
- fPidHF->SetCompat(kTRUE);\r
- }\r
-\r
-\r
- if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded\r
- isD0D0barPID[0]=0;\r
- isD0D0barPID[1]=0;\r
- }\r
- else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){\r
- if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded\r
- else isD0D0barPID[0]=0;// if K+ D0 excluded\r
- }\r
- /* else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){\r
- isD0D0barPID[0]=0;\r
- isD0D0barPID[1]=0;\r
- }\r
- */\r
- else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){ \r
- if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded\r
- else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded\r
- }\r
- else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){\r
- if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded\r
- else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded\r
- }\r
-\r
- if(fLowPt && d->Pt()<fPtLowPID){\r
- Double_t sigmaTPC[3]={3.,2.,0.};\r
- fPidHF->SetSigmaForTPC(sigmaTPC);\r
- // identify kaon\r
- combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);\r
-\r
- Double_t ptProng=aodtrack->P();\r
-\r
- if(ptProng<0.6){\r
- fPidHF->SetCompat(kFALSE);\r
- combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);\r
- fPidHF->SetCompat(kTRUE);\r
- }\r
-\r
- if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {\r
- combinedPID[daught][1]=0;\r
- }else{\r
- fPidHF->SetTOF(kFALSE);\r
- Double_t sigmaTPCpi[3]={3.,3.,0.};\r
- fPidHF->SetSigmaForTPC(sigmaTPCpi);\r
- combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);\r
- fPidHF->SetTOF(kTRUE);\r
- if(ptProng<0.8){\r
- Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");\r
- if(isTPCpion){\r
- combinedPID[daught][1]=1;\r
- }else{\r
- combinedPID[daught][1]=-1;\r
- }\r
- }\r
- }\r
-\r
- }\r
- fPidHF->SetSigmaForTPC(sigma_tmp);\r
- }// END OF LOOP ON DAUGHTERS\r
-\r
- if(!checkPIDInfo[0] && !checkPIDInfo[1]) {\r
- if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);\r
- return 0;\r
- }\r
-\r
-\r
- // FURTHER PID REQUEST (both daughter info is needed)\r
- if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){\r
- fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found\r
- if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);\r
- return 0;\r
- }\r
-\r
- if(fLowPt && d->Pt()<fPtLowPID){ \r
- if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){\r
- fWhyRejection=32;// reject cases where the Kaon is not identified\r
- fPidHF->SetSigmaForTPC(sigma_tmp);\r
- return 0;\r
- }\r
- }\r
- if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);\r
-\r
- // cout<<"Why? "<<fWhyRejection<<endl; \r
- return isD0D0barPID[0]+isD0D0barPID[1];\r
-}\r
-//---------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d) \r
-{\r
- // ############################################################\r
- //\r
- // Apply PID selection\r
- //\r
- //\r
- // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)\r
- //\r
- // d must be a AliAODRecoDecayHF2Prong object\r
- // returns 0 if both D0 and D0bar are rejectecd\r
- // 1 if D0 is accepted while D0bar is rejected\r
- // 2 if D0bar is accepted while D0 is rejected\r
- // 3 if both are accepted\r
- // fWhyRejection variable (not returned for the moment, print it if needed)\r
- // keeps some information on why a candidate has been \r
- // rejected according to the following (unfriendly?) scheme \r
- // if more rejection cases are considered interesting, just add numbers\r
- //\r
- // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant) \r
- // from 20 to 30: "detector" selection (PID acceptance) \r
- // 26: TPC refit\r
- // 27: ITS refit\r
- // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)\r
- //\r
- // from 30 to 40: PID selection\r
- // 31: no Kaon compatible tracks found between daughters\r
- // 32: no Kaon identified tracks found (strong sel. at low momenta)\r
- // 33: both mass hypotheses are rejected \r
- // \r
- // ############################################################\r
-\r
- if(!fUsePID) return 3;\r
- fWhyRejection=0;\r
- Int_t isD0D0barPID[2]={1,2};\r
- Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid\r
- Double_t tofSig,times[5];// used fot TOF pid\r
- Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters\r
- Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];\r
- Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value; \r
- // same for prong 2\r
- // values convention -1 = discarded \r
- // 0 = not identified (but compatible) || No PID (->hasPID flag)\r
- // 1 = identified\r
- // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both \r
- // Initial hypothesis: unknwon (but compatible) \r
- isKaonPionTOF[0][0]=0;\r
- isKaonPionTOF[0][1]=0;\r
- isKaonPionTOF[1][0]=0;\r
- isKaonPionTOF[1][1]=0;\r
- \r
- isKaonPionTPC[0][0]=0;\r
- isKaonPionTPC[0][1]=0;\r
- isKaonPionTPC[1][0]=0;\r
- isKaonPionTPC[1][1]=0;\r
- \r
- combinedPID[0][0]=0;\r
- combinedPID[0][1]=0;\r
- combinedPID[1][0]=0;\r
- combinedPID[1][1]=0;\r
- \r
- \r
- \r
- for(Int_t daught=0;daught<2;daught++){\r
- //Loop con prongs\r
- \r
- // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################\r
-\r
- AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught); \r
- \r
- if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){\r
- fWhyRejection=26;\r
- return 0;\r
- } \r
- if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){\r
- fWhyRejection=27;\r
- return 0;\r
- } \r
- \r
- AliAODPid *pid=aodtrack->GetDetPid();\r
- if(!pid) {\r
- //delete esdtrack;\r
- hasPID[daught]--;\r
- continue;\r
- }\r
- \r
- // ########### Step 1- Check of TPC and TOF response ####################\r
-\r
- Double_t ptrack=aodtrack->P();\r
- //#################### TPC PID #######################\r
- if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){\r
- // NO TPC PID INFO FOR THIS TRACK \r
- hasPID[daught]--;\r
- }\r
- else {\r
- static AliTPCPIDResponse theTPCpid;\r
- AliAODPid *pidObj = aodtrack->GetDetPid();\r
- Double_t ptProng=pidObj->GetTPCmomentum();\r
- nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);\r
- nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);\r
- //if(ptrack<0.6){\r
- if(ptProng<0.6){\r
- if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;\r
- else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;\r
- else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- }\r
- //else if(ptrack<.8){\r
- else if(ptProng<.8){\r
- if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;\r
- else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;\r
- else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- } \r
- else {\r
- // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;\r
- if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;\r
- if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- }\r
- }\r
- \r
- \r
- // ##### TOF PID: do not ask nothing for pion/protons ############\r
- if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){\r
- // NO TOF PID INFO FOR THIS TRACK \r
- hasPID[daught]--;\r
- }\r
- else{\r
- tofSig=pid->GetTOFsignal(); \r
- pid->GetIntegratedTimes(times);\r
- if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION\r
- if(TMath::Abs(tofSig-times[3])>3.*160.){\r
- isKaonPionTOF[daught][0]=-1;\r
- }\r
- else { \r
- if(ptrack<1.5){\r
- isKaonPionTOF[daught][0]=1;\r
- }\r
- }\r
- }\r
- \r
- //######### Step 2: COMBINE TOF and TPC PID ###############\r
- // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown\r
- combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];\r
- combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];\r
- \r
- \r
- //######### Step 3: USE PID INFO \r
- \r
- if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded\r
- isD0D0barPID[0]=0;\r
- isD0D0barPID[1]=0;\r
- }\r
- else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K\r
- if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded\r
- else isD0D0barPID[0]=0;// if K+ D0 excluded\r
- }\r
- else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject\r
- isD0D0barPID[0]=0;\r
- isD0D0barPID[1]=0;\r
- }\r
- else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){\r
- if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded\r
- else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded\r
- }\r
- else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){\r
- if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded\r
- else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded\r
- }\r
- \r
- // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################\r
- // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################\r
- // ############### NOT OPTIMIZED YET ###################################\r
- if(d->Pt()<2.){\r
- isKaonPionTPC[daught][0]=0;\r
- isKaonPionTPC[daught][1]=0;\r
- AliAODPid *pidObj = aodtrack->GetDetPid();\r
- Double_t ptProng=pidObj->GetTPCmomentum();\r
- //if(ptrack<0.6){\r
- if(ptProng<0.6){\r
- if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;\r
- else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;\r
- else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- }\r
- //else if(ptrack<.8){\r
- else if(ptProng<.8){\r
- if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;\r
- else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;\r
- else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- } \r
- else {\r
- if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;\r
- if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;\r
- }\r
- }\r
- \r
- }// END OF LOOP ON DAUGHTERS\r
- \r
- // FURTHER PID REQUEST (both daughter info is needed)\r
- if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){\r
- fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found\r
- return 0;\r
- }\r
- else if(hasPID[0]==0&&hasPID[1]==0){\r
- fWhyRejection=28;// reject cases in which no PID info is available \r
- return 0;\r
- }\r
- if(d->Pt()<2.){\r
- // request of K identification at low D0 pt\r
- combinedPID[0][0]=0;\r
- combinedPID[0][1]=0;\r
- combinedPID[1][0]=0;\r
- combinedPID[1][1]=0;\r
- \r
- combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];\r
- combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];\r
- combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];\r
- combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];\r
- \r
- if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){\r
- fWhyRejection=32;// reject cases where the Kaon is not identified\r
- return 0;\r
- }\r
- }\r
-\r
- // cout<<"Why? "<<fWhyRejection<<endl; \r
- return isD0D0barPID[0]+isD0D0barPID[1];\r
-}\r
-\r
-\r
-\r
-//---------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,\r
- Int_t selectionvalCand,\r
- Int_t selectionvalPID) const\r
-{\r
- //\r
- // This method combines the tracks, PID and cuts selection results\r
- //\r
- if(selectionvalTrack==0) return 0;\r
-\r
- Int_t returnvalue;\r
-\r
- switch(selectionvalPID) {\r
- case 0:\r
- returnvalue=0;\r
- break;\r
- case 1:\r
- returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);\r
- break;\r
- case 2:\r
- returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);\r
- break;\r
- case 3:\r
- returnvalue=selectionvalCand;\r
- break;\r
- default:\r
- returnvalue=0;\r
- break;\r
- }\r
-\r
- return returnvalue;\r
-}\r
-//----------------------------------------------------------------------------\r
-Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const\r
-{\r
- //\r
- // Note: this method is temporary\r
- // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters\r
- //\r
-\r
- //apply cuts\r
-\r
- Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;\r
- // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066\r
- // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution), \r
-\r
- Int_t returnvalue=3; //cut passed\r
- for(Int_t i=0;i<2/*prongs*/;i++){\r
- if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed\r
- }\r
- if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed\r
- if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed\r
- \r
- return returnvalue;\r
-}\r
-\r
-//----------------------------------------------\r
-void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)\r
-{\r
- //switch on candidate selection via AliKFparticle\r
- if(!useKF) return;\r
- if(useKF){\r
- fUseKF=useKF;\r
- Int_t nvarsKF=11;\r
- SetNVars(nvarsKF);\r
- TString varNamesKF[11]={"inv. mass [GeV]", \r
- "dca [cm]",\r
- "cosThetaStar", \r
- "pTK [GeV/c]",\r
- "pTPi [GeV/c]",\r
- "d0K [cm]",\r
- "d0Pi [cm]",\r
- "d0d0 [cm^2]",\r
- "cosThetaPoint"\r
- "DecayLength[cm]",\r
- "RedChi2"};\r
- Bool_t isUpperCutKF[11]={kTRUE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE,\r
- kTRUE};\r
- SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);\r
- SetGlobalIndex();\r
- Bool_t forOpt[11]={kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE,\r
- kFALSE,\r
- kFALSE,\r
- kTRUE,\r
- kTRUE,\r
- kFALSE,\r
- kFALSE};\r
- SetVarsForOpt(4,forOpt);\r
- }\r
- return;\r
-}\r
-\r
-\r
-void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {\r
- //\r
- //STANDARD CUTS USED FOR 2010 pp analysis \r
- //dca cut will be enlarged soon to 400 micron\r
- //\r
- \r
- SetName("D0toKpiCutsStandard");\r
- SetTitle("Standard Cuts for D0 analysis");\r
- \r
- // PILE UP REJECTION\r
- SetOptPileup(AliRDHFCuts::kRejectPileupEvent);\r
-\r
- // EVENT CUTS\r
- SetMinVtxContr(1);\r
- // MAX Z-VERTEX CUT\r
- SetMaxVtxZ(10.);\r
- \r
- // TRACKS ON SINGLE TRACKS\r
- AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");\r
- esdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
- esdTrackCuts->SetRequireTPCRefit(kTRUE);\r
- esdTrackCuts->SetRequireITSRefit(kTRUE);\r
- // esdTrackCuts->SetMinNClustersITS(4);\r
- esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
- esdTrackCuts->SetMinDCAToVertexXY(0.);\r
- esdTrackCuts->SetEtaRange(-0.8,0.8);\r
- esdTrackCuts->SetPtRange(0.3,1.e10);\r
- \r
- AddTrackCuts(esdTrackCuts);\r
-\r
- \r
- const Int_t nptbins =14;\r
- const Double_t ptmax = 9999.;\r
- const Int_t nvars=11;\r
- Float_t ptbins[nptbins+1];\r
- ptbins[0]=0.;\r
- ptbins[1]=0.5; \r
- ptbins[2]=1.;\r
- ptbins[3]=2.;\r
- ptbins[4]=3.;\r
- ptbins[5]=4.;\r
- ptbins[6]=5.;\r
- ptbins[7]=6.;\r
- ptbins[8]=7.;\r
- ptbins[9]=8.;\r
- ptbins[10]=12.;\r
- ptbins[11]=16.;\r
- ptbins[12]=20.;\r
- ptbins[13]=24.;\r
- ptbins[14]=ptmax;\r
-\r
- SetGlobalIndex(nvars,nptbins);\r
- SetPtBins(nptbins+1,ptbins);\r
- \r
- Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* pt<0.5*/\r
- {0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* 0.5<pt<1*/\r
- {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.80,0.,0.},/* 1<pt<2 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,0.},/* 2<pt<3 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 3<pt<4 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 4<pt<5 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 5<pt<6 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-7000.*1E-8,0.85,0.,0.},/* 7<pt<8 */\r
- {0.400,300.*1E-4,0.9,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.,0.},/* 8<pt<12 */\r
- {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,10000.*1E-8,0.85,0.,0.},/* 12<pt<16 */\r
- {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 16<pt<20 */\r
- {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 20<pt<24 */\r
- {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.}};/* pt>24 */\r
- \r
- \r
- //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts\r
- Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];\r
- for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];\r
- \r
- for (Int_t ibin=0;ibin<nptbins;ibin++){\r
- for (Int_t ivar = 0; ivar<nvars; ivar++){\r
- cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar]; \r
- }\r
- }\r
- \r
- SetCuts(nvars,nptbins,cutsMatrixTransposeStand);\r
- SetUseSpecialCuts(kTRUE);\r
- SetRemoveDaughtersFromPrim(kTRUE);\r
- \r
- for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];\r
- delete [] cutsMatrixTransposeStand;\r
- cutsMatrixTransposeStand=NULL;\r
-\r
- // PID SETTINGS\r
- AliAODPidHF* pidObj=new AliAODPidHF();\r
- //pidObj->SetName("pid4D0");\r
- Int_t mode=1;\r
- const Int_t nlims=2;\r
- Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]\r
- Bool_t compat=kTRUE; //effective only for this mode\r
- Bool_t asym=kTRUE;\r
- Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella\r
- pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC\r
- pidObj->SetMatch(mode);\r
- pidObj->SetPLimit(plims,nlims);\r
- pidObj->SetSigma(sigmas);\r
- pidObj->SetCompat(compat);\r
- pidObj->SetTPC(kTRUE);\r
- pidObj->SetTOF(kTRUE);\r
- pidObj->SetPCompatTOF(1.5);\r
- pidObj->SetSigmaForTPCCompat(3.);\r
- pidObj->SetSigmaForTOFCompat(3.);\r
- \r
-\r
- SetPidHF(pidObj);\r
- SetUsePID(kTRUE);\r
- SetUseDefaultPID(kFALSE);\r
-\r
- SetLowPt(kFALSE);\r
-\r
- PrintAll();\r
-\r
- delete pidObj;\r
- pidObj=NULL;\r
-\r
- return;\r
-\r
-}\r
-\r
-\r
-void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {\r
- //\r
- // Final CUTS USED FOR 2010 PbPb analysis of central events\r
- // REMEMBER TO EVENTUALLY SWITCH ON \r
- // SetUseAOD049(kTRUE);\r
- // \r
- \r
- SetName("D0toKpiCutsStandard");\r
- SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");\r
- \r
- \r
- // CENTRALITY SELECTION\r
- SetMinCentrality(0.);\r
- SetMaxCentrality(80.);\r
- SetUseCentrality(AliRDHFCuts::kCentV0M);\r
-\r
-\r
- \r
- // EVENT CUTS\r
- SetMinVtxContr(1);\r
- // MAX Z-VERTEX CUT\r
- SetMaxVtxZ(10.);\r
- // PILE UP\r
- SetOptPileup(AliRDHFCuts::kNoPileupSelection);\r
- // PILE UP REJECTION\r
- //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);\r
-\r
- // TRACKS ON SINGLE TRACKS\r
- AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");\r
- esdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
- esdTrackCuts->SetRequireTPCRefit(kTRUE);\r
- esdTrackCuts->SetRequireITSRefit(kTRUE);\r
- // esdTrackCuts->SetMinNClustersITS(4);\r
- esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
- esdTrackCuts->SetMinDCAToVertexXY(0.);\r
- esdTrackCuts->SetEtaRange(-0.8,0.8);\r
- esdTrackCuts->SetPtRange(0.7,1.e10);\r
-\r
- esdTrackCuts->SetMaxDCAToVertexXY(1.); \r
- esdTrackCuts->SetMaxDCAToVertexZ(1.);\r
- esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))"); \r
-\r
-\r
- AddTrackCuts(esdTrackCuts);\r
- // SPD k FIRST for Pt<3 GeV/c\r
- SetSelectCandTrackSPDFirst(kTRUE, 3); \r
-\r
- // CANDIDATE CUTS \r
- const Int_t nptbins =13;\r
- const Double_t ptmax = 9999.;\r
- const Int_t nvars=11;\r
- Float_t ptbins[nptbins+1];\r
- ptbins[0]=0.;\r
- ptbins[1]=0.5; \r
- ptbins[2]=1.;\r
- ptbins[3]=2.;\r
- ptbins[4]=3.;\r
- ptbins[5]=4.;\r
- ptbins[6]=5.;\r
- ptbins[7]=6.;\r
- ptbins[8]=8.;\r
- ptbins[9]=12.;\r
- ptbins[10]=16.;\r
- ptbins[11]=20.;\r
- ptbins[12]=24.;\r
- ptbins[13]=ptmax;\r
-\r
- SetGlobalIndex(nvars,nptbins);\r
- SetPtBins(nptbins+1,ptbins);\r
- SetMinPtCandidate(2.);\r
-\r
- Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.99,2.},/* pt<0.5*/\r
- {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.99,2.},/* 0.5<pt<1*/\r
- {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.995,2.},/* 1<pt<2 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-45000.*1E-8,0.95,0.998,7.},/* 2<pt<3 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */\r
- {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */\r
- {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,8.},/* 12<pt<16 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.995,6.},/* 16<pt<20 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.995,6.},/* 20<pt<24 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.99,2.}};/* pt>24 */\r
- \r
- \r
- //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts\r
- Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];\r
- for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];\r
- \r
- for (Int_t ibin=0;ibin<nptbins;ibin++){\r
- for (Int_t ivar = 0; ivar<nvars; ivar++){\r
- cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar]; \r
- }\r
- }\r
- \r
- SetCuts(nvars,nptbins,cutsMatrixTransposeStand);\r
- SetUseSpecialCuts(kTRUE);\r
- SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb\r
- for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];\r
- delete [] cutsMatrixTransposeStand;\r
- cutsMatrixTransposeStand=NULL;\r
- \r
- // PID SETTINGS\r
- AliAODPidHF* pidObj=new AliAODPidHF();\r
- //pidObj->SetName("pid4D0");\r
- Int_t mode=1;\r
- const Int_t nlims=2;\r
- Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]\r
- Bool_t compat=kTRUE; //effective only for this mode\r
- Bool_t asym=kTRUE;\r
- Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella\r
- pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC\r
- pidObj->SetMatch(mode);\r
- pidObj->SetPLimit(plims,nlims);\r
- pidObj->SetSigma(sigmas);\r
- pidObj->SetCompat(compat);\r
- pidObj->SetTPC(kTRUE);\r
- pidObj->SetTOF(kTRUE);\r
- pidObj->SetPCompatTOF(2.);\r
- pidObj->SetSigmaForTPCCompat(3.);\r
- pidObj->SetSigmaForTOFCompat(3.); \r
-\r
-\r
- SetPidHF(pidObj);\r
- SetUsePID(kTRUE);\r
- SetUseDefaultPID(kFALSE);\r
-\r
-\r
- PrintAll();\r
-\r
-\r
- delete pidObj;\r
- pidObj=NULL;\r
-\r
- return;\r
-\r
-}\r
-\r
-void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {\r
- // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%\r
-\r
- \r
- SetName("D0toKpiCutsStandard");\r
- SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");\r
- \r
- \r
- // CENTRALITY SELECTION\r
- SetMinCentrality(40.);\r
- SetMaxCentrality(80.);\r
- SetUseCentrality(AliRDHFCuts::kCentV0M);\r
- \r
- // EVENT CUTS\r
- SetMinVtxContr(1);\r
- // MAX Z-VERTEX CUT\r
- SetMaxVtxZ(10.);\r
- // PILE UP\r
- SetOptPileup(AliRDHFCuts::kNoPileupSelection);\r
- // PILE UP REJECTION\r
- //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);\r
-\r
- // TRACKS ON SINGLE TRACKS\r
- AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");\r
- esdTrackCuts->SetRequireSigmaToVertex(kFALSE);\r
- esdTrackCuts->SetRequireTPCRefit(kTRUE);\r
- esdTrackCuts->SetRequireITSRefit(kTRUE);\r
- // esdTrackCuts->SetMinNClustersITS(4);\r
- esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);\r
- esdTrackCuts->SetMinDCAToVertexXY(0.);\r
- esdTrackCuts->SetEtaRange(-0.8,0.8);\r
- esdTrackCuts->SetPtRange(0.5,1.e10);\r
-\r
- esdTrackCuts->SetMaxDCAToVertexXY(1.); \r
- esdTrackCuts->SetMaxDCAToVertexZ(1.);\r
- esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))"); \r
-\r
-\r
- AddTrackCuts(esdTrackCuts);\r
- // SPD k FIRST for Pt<3 GeV/c\r
- SetSelectCandTrackSPDFirst(kTRUE, 3); \r
-\r
- // CANDIDATE CUTS \r
- const Int_t nptbins =13;\r
- const Double_t ptmax = 9999.;\r
- const Int_t nvars=11;\r
- Float_t ptbins[nptbins+1];\r
- ptbins[0]=0.;\r
- ptbins[1]=0.5; \r
- ptbins[2]=1.;\r
- ptbins[3]=2.;\r
- ptbins[4]=3.;\r
- ptbins[5]=4.;\r
- ptbins[6]=5.;\r
- ptbins[7]=6.;\r
- ptbins[8]=8.;\r
- ptbins[9]=12.;\r
- ptbins[10]=16.;\r
- ptbins[11]=20.;\r
- ptbins[12]=24.;\r
- ptbins[13]=ptmax;\r
-\r
- SetGlobalIndex(nvars,nptbins);\r
- SetPtBins(nptbins+1,ptbins);\r
- SetMinPtCandidate(0.);\r
-\r
- Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.7,0.993,2.},/* pt<0.5*/\r
- {0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.85,0.993,2.},/* 0.5<pt<1*/\r
- {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.85,0.995,6.},/* 1<pt<2 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.991,5.},/* 2<pt<3 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.993,5.},/* 3<pt<4 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.995,5.},/* 4<pt<5 */\r
- {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */\r
- {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */\r
- {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.995,5.},/* 8<pt<12 */\r
- {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.83,0.995,4.},/* 12<pt<16 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.82,0.995,4.},/* 16<pt<20 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.81,0.995,4.},/* 20<pt<24 */\r
- {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.8,0.995,4.}};/* pt>24 */\r
- \r
- \r
- //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts\r
- Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];\r
- for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];\r
- \r
- for (Int_t ibin=0;ibin<nptbins;ibin++){\r
- for (Int_t ivar = 0; ivar<nvars; ivar++){\r
- cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar]; \r
- }\r
- }\r
- \r
- SetCuts(nvars,nptbins,cutsMatrixTransposeStand);\r
- SetUseSpecialCuts(kTRUE);\r
- SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb\r
- for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];\r
- delete [] cutsMatrixTransposeStand;\r
- cutsMatrixTransposeStand=NULL;\r
- \r
- // PID SETTINGS\r
- AliAODPidHF* pidObj=new AliAODPidHF();\r
- //pidObj->SetName("pid4D0");\r
- Int_t mode=1;\r
- const Int_t nlims=2;\r
- Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]\r
- Bool_t compat=kTRUE; //effective only for this mode\r
- Bool_t asym=kTRUE;\r
- Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella\r
- pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC\r
- pidObj->SetMatch(mode);\r
- pidObj->SetPLimit(plims,nlims);\r
- pidObj->SetSigma(sigmas);\r
- pidObj->SetCompat(compat);\r
- pidObj->SetTPC(kTRUE);\r
- pidObj->SetTOF(kTRUE);\r
- pidObj->SetPCompatTOF(2.);\r
- pidObj->SetSigmaForTPCCompat(3.);\r
- pidObj->SetSigmaForTOFCompat(3.); \r
-\r
- SetPidHF(pidObj);\r
- SetUsePID(kTRUE);\r
- SetUseDefaultPID(kFALSE);\r
-\r
- SetLowPt(kTRUE,2.);\r
-\r
- PrintAll();\r
-\r
-\r
- delete pidObj;\r
- pidObj=NULL;\r
-\r
- return;\r
- \r
-}\r
-\r
-\r
-void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {\r
- \r
- // Default 2010 PbPb cut object\r
- SetStandardCutsPbPb2010();\r
-\r
- //\r
- // Enable all 2011 PbPb run triggers\r
- // \r
- SetTriggerClass("");\r
- ResetMaskAndEnableMBTrigger();\r
- EnableCentralTrigger();\r
- EnableSemiCentralTrigger();\r
-\r
-}\r
+/**************************************************************************
+ * 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 D0->Kpi
+//
+// Author: A.Dainese, andrea.dainese@pd.infn.it
+/////////////////////////////////////////////////////////////
+
+#include <TDatabasePDG.h>
+#include <Riostream.h>
+
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliAODRecoDecayHF2Prong.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliAODPid.h"
+#include "AliTPCPIDResponse.h"
+#include "AliAODVertex.h"
+#include "AliKFVertex.h"
+#include "AliKFParticle.h"
+
+using std::cout;
+using std::endl;
+
+ClassImp(AliRDHFCutsD0toKpi)
+
+//--------------------------------------------------------------------------
+AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
+ AliRDHFCuts(name),
+ fUseSpecialCuts(kFALSE),
+ fLowPt(kTRUE),
+ fDefaultPID(kFALSE),
+ fUseKF(kFALSE),
+ fPtLowPID(2.),
+ fPtMaxSpecialCuts(9999.),
+ fmaxPtrackForPID(9999),
+ fCombPID(kFALSE),
+ fnSpecies(AliPID::kSPECIES),
+ fWeightsPositive(0),
+ fWeightsNegative(0),
+ fProbThreshold(0.5),
+ fBayesianStrategy(0),
+ fBayesianCondition(0)
+{
+ //
+ // Default Constructor
+ //
+ Int_t nvars=11;
+ SetNVars(nvars);
+ TString varNames[11]={"inv. mass [GeV]",
+ "dca [cm]",
+ "cosThetaStar",
+ "pTK [GeV/c]",
+ "pTPi [GeV/c]",
+ "d0K [cm]",
+ "d0Pi [cm]",
+ "d0d0 [cm^2]",
+ "cosThetaPoint",
+ "|cosThetaPointXY|",
+ "NormDecayLenghtXY"};
+ Bool_t isUpperCut[11]={kTRUE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kTRUE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kFALSE};
+ SetVarNames(nvars,varNames,isUpperCut);
+ Bool_t forOpt[11]={kFALSE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kFALSE,
+ kFALSE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE};
+ SetVarsForOpt(4,forOpt);
+ Float_t limits[2]={0,999999999.};
+ SetPtBins(2,limits);
+
+ fWeightsNegative = new Double_t[AliPID::kSPECIES];
+ fWeightsPositive = new Double_t[AliPID::kSPECIES];
+
+ for (Int_t i = 0; i<AliPID::kSPECIES; i++) {
+ fWeightsPositive[i] = 0.;
+ fWeightsNegative[i] = 0.;
+ }
+
+}
+//--------------------------------------------------------------------------
+AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
+ AliRDHFCuts(source),
+ fUseSpecialCuts(source.fUseSpecialCuts),
+ fLowPt(source.fLowPt),
+ fDefaultPID(source.fDefaultPID),
+ fUseKF(source.fUseKF),
+ fPtLowPID(source.fPtLowPID),
+ fPtMaxSpecialCuts(source.fPtMaxSpecialCuts),
+ fmaxPtrackForPID(source.fmaxPtrackForPID),
+ fCombPID(source.fCombPID),
+ fnSpecies(source.fnSpecies),
+ fWeightsPositive(source.fWeightsPositive),
+ fWeightsNegative(source.fWeightsNegative),
+ fProbThreshold(source.fProbThreshold),
+ fBayesianStrategy(source.fBayesianStrategy),
+ fBayesianCondition(source.fBayesianCondition)
+{
+ //
+ // Copy constructor
+ //
+
+}
+//--------------------------------------------------------------------------
+AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
+{
+ //
+ // assignment operator
+ //
+ if(&source == this) return *this;
+
+ AliRDHFCuts::operator=(source);
+ fUseSpecialCuts=source.fUseSpecialCuts;
+ fLowPt=source.fLowPt;
+ fDefaultPID=source.fDefaultPID;
+ fUseKF=source.fUseKF;
+ fPtLowPID=source.fPtLowPID;
+ fPtMaxSpecialCuts=source.fPtMaxSpecialCuts;
+ fmaxPtrackForPID=source.fmaxPtrackForPID;
+ fCombPID = source.fCombPID;
+ fnSpecies = source.fnSpecies;
+ fWeightsPositive = source.fWeightsPositive;
+ fWeightsNegative = source.fWeightsNegative;
+ fProbThreshold = source.fProbThreshold;
+ fBayesianStrategy = source.fBayesianStrategy;
+ fBayesianCondition = source.fBayesianCondition;
+ return *this;
+}
+
+
+
+//---------------------------------------------------------------------------
+AliRDHFCutsD0toKpi::~AliRDHFCutsD0toKpi()
+{
+//
+// Destructor
+//
+ if (fWeightsNegative) {
+ delete[] fWeightsNegative;
+ fWeightsNegative = 0;
+ }
+ if (fWeightsPositive) {
+ delete[] fWeightsNegative;
+ fWeightsNegative = 0;
+ }
+}
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
+ //
+ // Fills in vars the values of the variables
+ //
+
+ if(nvars!=fnVarsForOpt) {
+ printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
+ return;
+ }
+
+ AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
+
+ //recalculate vertex w/o daughters
+ Bool_t cleanvtx=kFALSE;
+ AliAODVertex *origownvtx=0x0;
+ if(fRemoveDaughtersFromPrimary) {
+ if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
+ cleanvtx=kTRUE;
+ if(!RecalcOwnPrimaryVtx(dd,aod)) {
+ CleanOwnPrimaryVtx(dd,aod,origownvtx);
+ cleanvtx=kFALSE;
+ }
+ }
+
+ Int_t iter=-1;
+ if(fVarsForOpt[0]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==211) {
+ vars[iter]=dd->InvMassD0();
+ } else {
+ vars[iter]=dd->InvMassD0bar();
+ }
+ }
+ if(fVarsForOpt[1]){
+ iter++;
+ vars[iter]=dd->GetDCA();
+ }
+ if(fVarsForOpt[2]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==211) {
+ vars[iter] = dd->CosThetaStarD0();
+ } else {
+ vars[iter] = dd->CosThetaStarD0bar();
+ }
+ }
+ if(fVarsForOpt[3]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==321) {
+ vars[iter]=dd->PtProng(0);
+ }
+ else{
+ vars[iter]=dd->PtProng(1);
+ }
+ }
+ if(fVarsForOpt[4]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==211) {
+ vars[iter]=dd->PtProng(0);
+ }
+ else{
+ vars[iter]=dd->PtProng(1);
+ }
+ }
+ if(fVarsForOpt[5]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==321) {
+ vars[iter]=dd->Getd0Prong(0);
+ }
+ else{
+ vars[iter]=dd->Getd0Prong(1);
+ }
+ }
+ if(fVarsForOpt[6]){
+ iter++;
+ if(TMath::Abs(pdgdaughters[0])==211) {
+ vars[iter]=dd->Getd0Prong(0);
+ }
+ else{
+ vars[iter]=dd->Getd0Prong(1);
+ }
+ }
+ if(fVarsForOpt[7]){
+ iter++;
+ vars[iter]= dd->Prodd0d0();
+ }
+ if(fVarsForOpt[8]){
+ iter++;
+ vars[iter]=dd->CosPointingAngle();
+ }
+
+ if(fVarsForOpt[9]){
+ iter++;
+ vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
+ }
+
+ if(fVarsForOpt[10]){
+ iter++;
+ vars[iter]=dd->NormalizedDecayLengthXY();
+ }
+
+ if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
+
+ return;
+}
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
+ //
+ // Apply selection
+ //
+
+
+ fIsSelectedCuts=0;
+ fIsSelectedPID=0;
+
+ if(!fCutsRD){
+ cout<<"Cut matrice not inizialized. Exit..."<<endl;
+ return 0;
+ }
+ //PrintAll();
+ AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
+ if(!d){
+ cout<<"AliAODRecoDecayHF2Prong null"<<endl;
+ return 0;
+ }
+
+ if(fKeepSignalMC) if(IsSignalMC(d,aod,421)) return 3;
+
+ Double_t ptD=d->Pt();
+ if(ptD<fMinPtCand) return 0;
+ if(ptD>fMaxPtCand) return 0;
+
+ if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
+
+ // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
+ Int_t returnvaluePID=3;
+ Int_t returnvalueCuts=3;
+
+ // selection on candidate
+ if(selectionLevel==AliRDHFCuts::kAll ||
+ selectionLevel==AliRDHFCuts::kCandidate) {
+
+ if(!fUseKF) {
+
+ //recalculate vertex w/o daughters
+ AliAODVertex *origownvtx=0x0;
+ if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
+ if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+ if(!RecalcOwnPrimaryVtx(d,aod)) {
+ CleanOwnPrimaryVtx(d,aod,origownvtx);
+ return 0;
+ }
+ }
+
+ if(fUseMCVertex) {
+ if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+ if(!SetMCPrimaryVtx(d,aod)) {
+ CleanOwnPrimaryVtx(d,aod,origownvtx);
+ return 0;
+ }
+ }
+
+ Double_t pt=d->Pt();
+
+ Int_t okD0=0,okD0bar=0;
+
+ Int_t ptbin=PtBin(pt);
+ if (ptbin==-1) {
+ CleanOwnPrimaryVtx(d,aod,origownvtx);
+ return 0;
+ }
+
+ Double_t mD0,mD0bar,ctsD0,ctsD0bar;
+ okD0=1; okD0bar=1;
+
+ Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+
+ d->InvMassD0(mD0,mD0bar);
+ if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
+ if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+
+ if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
+ if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+
+ if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
+ TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
+ if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
+ TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+
+ d->CosThetaStarD0(ctsD0,ctsD0bar);
+ if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
+ if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+
+ if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ Double_t normalDecayLengXY=d->NormalizedDecayLengthXY();
+ if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ if (returnvalueCuts!=0) {
+ if (okD0) returnvalueCuts=1; //cuts passed as D0
+ if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
+ if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
+ }
+
+ // call special cuts
+ Int_t special=1;
+ if(fUseSpecialCuts && (pt<fPtMaxSpecialCuts)) special=IsSelectedSpecialCuts(d);
+ if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+ // unset recalculated primary vertex when not needed any more
+ CleanOwnPrimaryVtx(d,aod,origownvtx);
+
+ } else {
+ // go to selection with Kalman vertexing, if requested
+ returnvalueCuts = IsSelectedKF(d,aod);
+ }
+
+ fIsSelectedCuts=returnvalueCuts;
+ if(!returnvalueCuts) return 0;
+ }
+
+
+ // selection on PID
+ if(selectionLevel==AliRDHFCuts::kAll ||
+ selectionLevel==AliRDHFCuts::kCandidate ||
+ selectionLevel==AliRDHFCuts::kPID) {
+ if (!fCombPID) {
+ returnvaluePID = IsSelectedPID(d);
+ fIsSelectedPID=returnvaluePID;
+ if(!returnvaluePID) return 0;
+ } else {
+ returnvaluePID = IsSelectedCombPID(d);
+ if(!returnvaluePID) return 0;
+ }
+ }
+
+
+
+ Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
+
+ if(!returnvalueComb) return 0;
+
+ // selection on daughter tracks
+ if(selectionLevel==AliRDHFCuts::kAll ||
+ selectionLevel==AliRDHFCuts::kTracks) {
+ if(!AreDaughtersSelected(d)) return 0;
+ }
+
+ // cout<<"Pid = "<<returnvaluePID<<endl;
+ return returnvalueComb;
+}
+
+//------------------------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
+ AliAODEvent* aod) const {
+ //
+ // Apply selection using KF-vertexing
+ //
+
+ AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
+ AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
+
+ if(!track0 || !track1) {
+ cout<<"one or two D0 daughters missing!"<<endl;
+ return 0;
+ }
+
+ // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
+ Int_t returnvalueCuts=3;
+
+ // candidate selection with AliKF
+ AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
+
+ Int_t okD0=0,okD0bar=0;
+ okD0=1; okD0bar=1;
+
+ // convert tracks into AliKFParticles
+
+ AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
+ AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
+ AliKFParticle posPiKF(*track0,211); // pos pion kandidate
+ AliKFParticle posKKF(*track0,321); // pos kaon kandidate
+
+ // build D0 candidates
+
+ AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
+ AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
+
+ // create kf primary vertices
+
+ AliAODVertex *vtx1 = aod->GetPrimaryVertex();
+ AliKFVertex primVtx(*vtx1);
+ AliKFVertex aprimVtx(*vtx1);
+
+ if(primVtx.GetNContributors()<=0) okD0 = 0;
+ if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ // calculate mass
+
+ Double_t d0mass = d0c.GetMass();
+ Double_t ad0mass = ad0c.GetMass();
+
+ // calculate P of D0 and D0bar
+ Double_t d0P = d0c.GetP();
+ Double_t d0Px = d0c.GetPx();
+ Double_t d0Py = d0c.GetPy();
+ Double_t d0Pz = d0c.GetPz();
+ Double_t ad0P = ad0c.GetP();
+ Double_t ad0Px = ad0c.GetPx();
+ Double_t ad0Py = ad0c.GetPy();
+ Double_t ad0Pz = ad0c.GetPz();
+
+ //calculate Pt of D0 and D0bar
+
+ Double_t pt=d0c.GetPt();
+ Double_t apt=ad0c.GetPt();
+
+ // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
+
+ if(track0->GetUsedForPrimVtxFit()) {
+ primVtx -= posPiKF;
+ aprimVtx -= posKKF;
+ }
+
+ if(track1->GetUsedForPrimVtxFit()) {
+ primVtx -= negKKF;
+ aprimVtx -= negPiKF;
+ }
+
+ primVtx += d0c;
+ aprimVtx += ad0c;
+
+ if(primVtx.GetNContributors()<=0) okD0 = 0;
+ if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ //calculate cut variables
+
+ // calculate impact params of daughters w.r.t recalculated vertices
+
+ Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
+ Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
+ Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
+ Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
+
+ // calculate Product of Impact Params
+
+ Double_t prodParam = impactPi*impactKa;
+ Double_t aprodParam = aimpactPi*aimpactKa;
+
+ // calculate cosine of pointing angles
+
+ TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
+ TVector3 fline(d0c.GetX()-primVtx.GetX(),
+ d0c.GetY()-primVtx.GetY(),
+ d0c.GetZ()-primVtx.GetZ());
+ Double_t pta = mom.Angle(fline);
+ Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
+
+ TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
+ TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
+ ad0c.GetY()-aprimVtx.GetY(),
+ ad0c.GetZ()-aprimVtx.GetZ());
+ Double_t apta = amom.Angle(afline);
+ Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
+
+ // calculate P of Pions at Decay Position of D0 and D0bar candidates
+ negKKF.TransportToParticle(d0c);
+ posPiKF.TransportToParticle(d0c);
+ posKKF.TransportToParticle(ad0c);
+ negPiKF.TransportToParticle(ad0c);
+
+ Double_t pxPi = posPiKF.GetPx();
+ Double_t pyPi = posPiKF.GetPy();
+ Double_t pzPi = posPiKF.GetPz();
+ Double_t ptPi = posPiKF.GetPt();
+
+ Double_t apxPi = negPiKF.GetPx();
+ Double_t apyPi = negPiKF.GetPy();
+ Double_t apzPi = negPiKF.GetPz();
+ Double_t aptPi = negPiKF.GetPt();
+
+ // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
+
+ Double_t ptK = negKKF.GetPt();
+ Double_t aptK = posKKF.GetPt();
+
+ //calculate cos(thetastar)
+ Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+ Double_t massp[2];
+ massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
+ massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+ Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
+ -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
+
+ // cos(thetastar) for D0 and Pion
+
+ Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
+ Double_t beta = d0P/d0E;
+ Double_t gamma = d0E/massvtx;
+ TVector3 momPi(pxPi,pyPi,pzPi);
+ TVector3 momTot(d0Px,d0Py,d0Pz);
+ Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
+ Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
+
+ // cos(thetastar) for D0bar and Pion
+
+ Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
+ Double_t abeta = ad0P/ad0E;
+ Double_t agamma = ad0E/massvtx;
+ TVector3 amomPi(apxPi,apyPi,apzPi);
+ TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
+ Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
+ Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
+
+ // calculate reduced Chi2 for the full D0 fit
+ d0c.SetProductionVertex(primVtx);
+ ad0c.SetProductionVertex(aprimVtx);
+ negKKF.SetProductionVertex(d0c);
+ posPiKF.SetProductionVertex(d0c);
+ posKKF.SetProductionVertex(ad0c);
+ negPiKF.SetProductionVertex(ad0c);
+ d0c.TransportToProductionVertex();
+ ad0c.TransportToProductionVertex();
+
+ // calculate the decay length
+ Double_t decayLengthD0 = d0c.GetDecayLength();
+ Double_t adecayLengthD0 = ad0c.GetDecayLength();
+
+ Double_t chi2D0 = 50.;
+ if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
+ chi2D0 = d0c.GetChi2()/d0c.GetNDF();
+ }
+
+ Double_t achi2D0 = 50.;
+ if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
+ achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
+ }
+
+ // Get the Pt-bins
+ Int_t ptbin=PtBin(pt);
+ Int_t aptbin=PtBin(apt);
+
+ if(ptbin < 0) okD0 = 0;
+ if(aptbin < 0) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
+ if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+
+ if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
+ TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
+
+ if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
+ TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
+
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ // for the moment via the standard method due to bug in AliKF
+ if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
+ if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+
+ if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
+ if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+
+ if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
+ if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(prodParam > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
+ if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(cosP < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
+ if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(chi2D0 > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
+ if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(decayLengthD0 < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
+ if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
+ if(!okD0 && !okD0bar) returnvalueCuts=0;
+
+ if(returnvalueCuts!=0) {
+ if(okD0) returnvalueCuts=1; //cuts passed as D0
+ if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
+ if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
+ }
+
+ return returnvalueCuts;
+}
+
+//---------------------------------------------------------------------------
+
+Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
+{
+ //
+ // Checking if D0 is in fiducial acceptance region
+ //
+
+ if(fMaxRapidityCand>-998.){
+ if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
+ else return kTRUE;
+ }
+
+ if(pt > 5.) {
+ // applying cut for pt > 5 GeV
+ AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",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 D0 = %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 AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
+{
+ // ############################################################
+ //
+ // Apply PID selection
+ //
+ //
+ // ############################################################
+
+ if(!fUsePID) return 3;
+ if(fDefaultPID) return IsSelectedPIDdefault(d);
+ if (fCombPID) return IsSelectedCombPID(d); //to use Bayesian method if applicable
+ fWhyRejection=0;
+ Int_t isD0D0barPID[2]={1,2};
+ Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
+ // same for prong 2
+ // values convention -1 = discarded
+ // 0 = not identified (but compatible) || No PID (->hasPID flag)
+ // 1 = identified
+ // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
+ // Initial hypothesis: unknwon (but compatible)
+ combinedPID[0][0]=0; // prima figlia, Kaon
+ combinedPID[0][1]=0; // prima figlia, pione
+ combinedPID[1][0]=0; // seconda figlia, Kaon
+ combinedPID[1][1]=0; // seconda figlia, pion
+
+ Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
+ Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
+ Bool_t isTOFused=fPidHF->GetTOF(),isCompat=fPidHF->GetCompat();
+ for(Int_t daught=0;daught<2;daught++){
+ //Loop con prongs
+ AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
+ if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
+
+ if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
+ checkPIDInfo[daught]=kFALSE;
+ continue;
+ }
+ Double_t pProng=aodtrack->P();
+
+ // identify kaon
+ if(pProng<fmaxPtrackForPID){
+ combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
+ }
+ // identify pion
+ if(pProng<fmaxPtrackForPID){
+ if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
+ combinedPID[daught][1]=0;
+ }else{
+ fPidHF->SetTOF(kFALSE);
+ combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
+ if(isTOFused)fPidHF->SetTOF(kTRUE);
+ if(isCompat)fPidHF->SetCompat(kTRUE);
+ }
+ }
+
+ if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
+ else isD0D0barPID[0]=0;// if K+ D0 excluded
+ }
+ /* else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ */
+ else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
+ else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
+ }
+ else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
+ else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
+ }
+
+ if(fLowPt && d->Pt()<fPtLowPID){
+ Double_t sigmaTPC[3]={3.,2.,0.};
+ fPidHF->SetSigmaForTPC(sigmaTPC);
+ // identify kaon
+ combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
+
+ if(pProng<0.6){
+ fPidHF->SetCompat(kFALSE);
+ combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
+ fPidHF->SetCompat(kTRUE);
+ }
+
+ if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
+ combinedPID[daught][1]=0;
+ }else{
+ fPidHF->SetTOF(kFALSE);
+ Double_t sigmaTPCpi[3]={3.,3.,0.};
+ fPidHF->SetSigmaForTPC(sigmaTPCpi);
+ combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
+ fPidHF->SetTOF(kTRUE);
+ if(pProng<0.8){
+ Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
+ if(isTPCpion){
+ combinedPID[daught][1]=1;
+ }else{
+ combinedPID[daught][1]=-1;
+ }
+ }
+ }
+
+ }
+ fPidHF->SetSigmaForTPC(sigma_tmp);
+ }// END OF LOOP ON DAUGHTERS
+
+ if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
+ if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
+ return 0;
+ }
+
+
+ // FURTHER PID REQUEST (both daughter info is needed)
+ if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
+ fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
+ if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
+ return 0;
+ }
+
+ if(fLowPt && d->Pt()<fPtLowPID){
+ if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
+ fWhyRejection=32;// reject cases where the Kaon is not identified
+ fPidHF->SetSigmaForTPC(sigma_tmp);
+ return 0;
+ }
+ }
+ if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
+
+ // cout<<"Why? "<<fWhyRejection<<endl;
+ return isD0D0barPID[0]+isD0D0barPID[1];
+}
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
+{
+ // ############################################################
+ //
+ // Apply PID selection
+ //
+ //
+ // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
+ //
+ // d must be a AliAODRecoDecayHF2Prong object
+ // returns 0 if both D0 and D0bar are rejectecd
+ // 1 if D0 is accepted while D0bar is rejected
+ // 2 if D0bar is accepted while D0 is rejected
+ // 3 if both are accepted
+ // fWhyRejection variable (not returned for the moment, print it if needed)
+ // keeps some information on why a candidate has been
+ // rejected according to the following (unfriendly?) scheme
+ // if more rejection cases are considered interesting, just add numbers
+ //
+ // TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
+ // from 20 to 30: "detector" selection (PID acceptance)
+ // 26: TPC refit
+ // 27: ITS refit
+ // 28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
+ //
+ // from 30 to 40: PID selection
+ // 31: no Kaon compatible tracks found between daughters
+ // 32: no Kaon identified tracks found (strong sel. at low momenta)
+ // 33: both mass hypotheses are rejected
+ //
+ // ############################################################
+
+ if(!fUsePID) return 3;
+ fWhyRejection=0;
+ Int_t isD0D0barPID[2]={1,2};
+ Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
+ Double_t tofSig,times[5];// used fot TOF pid
+ Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
+ Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
+ Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
+ // same for prong 2
+ // values convention -1 = discarded
+ // 0 = not identified (but compatible) || No PID (->hasPID flag)
+ // 1 = identified
+ // PID search: pion (TPC) or not K (TOF), Kaon hypothesis for both
+ // Initial hypothesis: unknwon (but compatible)
+ isKaonPionTOF[0][0]=0;
+ isKaonPionTOF[0][1]=0;
+ isKaonPionTOF[1][0]=0;
+ isKaonPionTOF[1][1]=0;
+
+ isKaonPionTPC[0][0]=0;
+ isKaonPionTPC[0][1]=0;
+ isKaonPionTPC[1][0]=0;
+ isKaonPionTPC[1][1]=0;
+
+ combinedPID[0][0]=0;
+ combinedPID[0][1]=0;
+ combinedPID[1][0]=0;
+ combinedPID[1][1]=0;
+
+
+
+ for(Int_t daught=0;daught<2;daught++){
+ //Loop con prongs
+
+ // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
+
+ AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
+
+ if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
+ fWhyRejection=26;
+ return 0;
+ }
+ if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
+ fWhyRejection=27;
+ return 0;
+ }
+
+ AliAODPid *pid=aodtrack->GetDetPid();
+ if(!pid) {
+ //delete esdtrack;
+ hasPID[daught]--;
+ continue;
+ }
+
+ // ########### Step 1- Check of TPC and TOF response ####################
+
+ Double_t ptrack=aodtrack->P();
+ //#################### TPC PID #######################
+ if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
+ // NO TPC PID INFO FOR THIS TRACK
+ hasPID[daught]--;
+ }
+ else {
+ static AliTPCPIDResponse theTPCpid;
+ AliAODPid *pidObj = aodtrack->GetDetPid();
+ Double_t ptProng=pidObj->GetTPCmomentum();
+ nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
+ nsigmaTPCK = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
+ //if(ptrack<0.6){
+ if(ptProng<0.6){
+ if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ //else if(ptrack<.8){
+ else if(ptProng<.8){
+ if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else {
+ // if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
+ if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
+ if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ }
+
+
+ // ##### TOF PID: do not ask nothing for pion/protons ############
+ if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
+ // NO TOF PID INFO FOR THIS TRACK
+ hasPID[daught]--;
+ }
+ else{
+ tofSig=pid->GetTOFsignal();
+ pid->GetIntegratedTimes(times);
+ if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
+ if(TMath::Abs(tofSig-times[3])>3.*160.){
+ isKaonPionTOF[daught][0]=-1;
+ }
+ else {
+ if(ptrack<1.5){
+ isKaonPionTOF[daught][0]=1;
+ }
+ }
+ }
+
+ //######### Step 2: COMBINE TOF and TPC PID ###############
+ // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
+ combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
+ combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
+
+
+ //######### Step 3: USE PID INFO
+
+ if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
+ else isD0D0barPID[0]=0;// if K+ D0 excluded
+ }
+ else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject
+ isD0D0barPID[0]=0;
+ isD0D0barPID[1]=0;
+ }
+ else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
+ else isD0D0barPID[0]=0;// not a D0 if K+ or if pi+ excluded
+ }
+ else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
+ if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
+ else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
+ }
+
+ // ########## ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification ###############################
+ // ########## more tolerant criteria for single particle ID-> more selective criteria for D0 ##############################
+ // ############### NOT OPTIMIZED YET ###################################
+ if(d->Pt()<2.){
+ isKaonPionTPC[daught][0]=0;
+ isKaonPionTPC[daught][1]=0;
+ AliAODPid *pidObj = aodtrack->GetDetPid();
+ Double_t ptProng=pidObj->GetTPCmomentum();
+ //if(ptrack<0.6){
+ if(ptProng<0.6){
+ if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ //else if(ptrack<.8){
+ else if(ptProng<.8){
+ if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
+ else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
+ else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ else {
+ if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
+ if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
+ }
+ }
+
+ }// END OF LOOP ON DAUGHTERS
+
+ // FURTHER PID REQUEST (both daughter info is needed)
+ if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
+ fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
+ return 0;
+ }
+ else if(hasPID[0]==0&&hasPID[1]==0){
+ fWhyRejection=28;// reject cases in which no PID info is available
+ return 0;
+ }
+ if(d->Pt()<2.){
+ // request of K identification at low D0 pt
+ combinedPID[0][0]=0;
+ combinedPID[0][1]=0;
+ combinedPID[1][0]=0;
+ combinedPID[1][1]=0;
+
+ combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
+ combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
+ combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
+ combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
+
+ if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
+ fWhyRejection=32;// reject cases where the Kaon is not identified
+ return 0;
+ }
+ }
+
+ // cout<<"Why? "<<fWhyRejection<<endl;
+ return isD0D0barPID[0]+isD0D0barPID[1];
+}
+
+
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
+ Int_t selectionvalCand,
+ Int_t selectionvalPID) const
+{
+ //
+ // This method combines the tracks, PID and cuts selection results
+ //
+ if(selectionvalTrack==0) return 0;
+
+ Int_t returnvalue;
+
+ switch(selectionvalPID) {
+ case 0:
+ returnvalue=0;
+ break;
+ case 1:
+ returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
+ break;
+ case 2:
+ returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
+ break;
+ case 3:
+ returnvalue=selectionvalCand;
+ break;
+ default:
+ returnvalue=0;
+ break;
+ }
+
+ return returnvalue;
+}
+//----------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
+{
+ //
+ // Note: this method is temporary
+ // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
+ //
+
+ //apply cuts
+
+ Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
+ // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
+ // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
+
+ Int_t returnvalue=3; //cut passed
+ for(Int_t i=0;i<2/*prongs*/;i++){
+ if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
+ }
+ if(d->DecayLength2()<decLengthCut*decLengthCut) return 0; //decLengthCut not passed
+ if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut) return 0; //decLengthCut not passed
+
+ return returnvalue;
+}
+
+//----------------------------------------------
+void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
+{
+ //switch on candidate selection via AliKFparticle
+ if(!useKF) return;
+ if(useKF){
+ fUseKF=useKF;
+ Int_t nvarsKF=11;
+ SetNVars(nvarsKF);
+ TString varNamesKF[11]={"inv. mass [GeV]",
+ "dca [cm]",
+ "cosThetaStar",
+ "pTK [GeV/c]",
+ "pTPi [GeV/c]",
+ "d0K [cm]",
+ "d0Pi [cm]",
+ "d0d0 [cm^2]",
+ "cosThetaPoint"
+ "DecayLength[cm]",
+ "RedChi2"};
+ Bool_t isUpperCutKF[11]={kTRUE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kTRUE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kTRUE};
+ SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
+ SetGlobalIndex();
+ Bool_t forOpt[11]={kFALSE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE,
+ kFALSE,
+ kFALSE,
+ kTRUE,
+ kTRUE,
+ kFALSE,
+ kFALSE};
+ SetVarsForOpt(4,forOpt);
+ }
+ return;
+}
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
+ //
+ //STANDARD CUTS USED FOR 2010 pp analysis
+ //dca cut will be enlarged soon to 400 micron
+ //
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis");
+
+ // PILE UP REJECTION
+ SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ // EVENT CUTS
+ SetMinVtxContr(1);
+ // MAX Z-VERTEX CUT
+ SetMaxVtxZ(10.);
+
+ // TRACKS ON SINGLE TRACKS
+ AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ // esdTrackCuts->SetMinNClustersITS(4);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+
+ const Int_t nptbins =14;
+ const Double_t ptmax = 9999.;
+ const Int_t nvars=11;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=7.;
+ ptbins[9]=8.;
+ ptbins[10]=12.;
+ ptbins[11]=16.;
+ ptbins[12]=20.;
+ ptbins[13]=24.;
+ ptbins[14]=ptmax;
+
+ SetGlobalIndex(nvars,nptbins);
+ SetPtBins(nptbins+1,ptbins);
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* pt<0.5*/
+ {0.400,350.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,0.},/* 0.5<pt<1*/
+ {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.80,0.,0.},/* 1<pt<2 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,0.},/* 2<pt<3 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 3<pt<4 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 4<pt<5 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 5<pt<6 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-7000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
+ {0.400,300.*1E-4,0.9,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.,0.},/* 8<pt<12 */
+ {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,10000.*1E-8,0.85,0.,0.},/* 12<pt<16 */
+ {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 16<pt<20 */
+ {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.},/* 20<pt<24 */
+ {0.400,300.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,999999.*1E-8,0.85,0.,0.}};/* pt>24 */
+
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ SetUseSpecialCuts(kTRUE);
+ SetRemoveDaughtersFromPrim(kTRUE);
+
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+ cutsMatrixTransposeStand=NULL;
+
+ // PID SETTINGS
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ //pidObj->SetName("pid4D0");
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetPCompatTOF(1.5);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+
+ SetPidHF(pidObj);
+ SetUsePID(kTRUE);
+ SetUseDefaultPID(kFALSE);
+
+ SetLowPt(kFALSE);
+
+ PrintAll();
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+
+}
+
+//___________________________________________________________________________
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2010vsMult() {
+ //
+ // Cuts for 2010 pp 7 TeV data analysis in Ntracklets bins (vs multiplicity)
+ //
+ SetName("D0toKpiCuts");
+ SetTitle("Cuts for D0 analysis in 2010-data pp 7 TeV vs multiplicity");
+
+ //
+ // Track cuts
+ //
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ //default
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+ AliESDtrackCuts::kAny);
+ // default is kBoth, otherwise kAny
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+
+
+ //
+ // Cut values per pt bin
+ //
+ const Int_t nvars=11;
+ const Int_t nptbins=14;
+ Float_t* ptbins;
+ ptbins=new Float_t[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=7.;
+ ptbins[9]=8.;
+ ptbins[10]=12.;
+ ptbins[11]=16.;
+ ptbins[12]=20.;
+ ptbins[13]=24.;
+ ptbins[14]=9999.;
+
+ SetPtBins(nptbins+1,ptbins);
+
+ //setting cut values
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* pt<0.5*/
+ {0.400,350.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.80,0.,3.},/* 0.5<pt<1*/
+ {0.400,300.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-30000.*1E-8,0.80,0.,4.},/* 1<pt<2 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.85,0.,4.},/* 2<pt<3 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.85,0.,4.},/* 3<pt<4 */
+ {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 4<pt<5 */
+ {0.400,300.*1E-4,0.75,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.875,0.,0.},/* 5<pt<6 */
+ {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<7 */
+ {0.400,400.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 7<pt<8 */
+ {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
+
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+
+ SetUseSpecialCuts(kTRUE);
+ SetMaximumPtSpecialCuts(8.);
+
+ //Do recalculate the vertex
+ SetRemoveDaughtersFromPrim(kTRUE);
+
+ //
+ // Pid settings
+ //
+ Bool_t pidflag=kTRUE;
+ SetUsePID(pidflag);
+ if(pidflag) cout<<"PID is used"<<endl;
+ else cout<<"PID is not used"<<endl;
+ //
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetPCompatTOF(1.5);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+ // pidObj->SetOldPid(kFALSE);
+ SetPidHF(pidObj);
+
+ SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
+
+ SetLowPt(kFALSE);
+
+ //activate pileup rejection (for pp)
+ SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+
+ PrintAll();
+
+ delete pidObj;
+ pidObj=NULL;
+ delete [] ptbins;
+ ptbins=NULL;
+
+ return;
+}
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
+ //
+ // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
+ //
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
+
+ //
+ // Track cuts
+ //
+ AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ //default
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
+ AliESDtrackCuts::kAny);
+ // default is kBoth, otherwise kAny
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetPtRange(0.3,1.e10);
+
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+
+
+ const Int_t nvars=11;
+ const Int_t nptbins=13;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=8.;
+ ptbins[9]=12.;
+ ptbins[10]=16.;
+ ptbins[11]=20.;
+ ptbins[12]=24.;
+ ptbins[13]=9999.;
+
+ SetPtBins(nptbins+1,ptbins);
+
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* pt<0.5*/
+ {0.400,0.04,0.75,0.3,0.3,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 0.5<pt<1*/
+ {0.400,0.03,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.8,0.,0.},/* 1<pt<2 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0003,0.9,0.,0.},/* 2<pt<3 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0002,0.9,0.,0.},/* 3<pt<4 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00015,0.9,0.,0.},/* 4<pt<5 */
+ {0.400,0.03,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.0001,0.9,0.,0.},/* 5<pt<6 */
+ {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
+ {0.400,0.06,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.00001,0.85,0.,0.},/* 8<pt<12 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
+ {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+
+
+ //pid settings
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetOldPid(kTRUE);
+ SetPidHF(pidObj);
+
+ SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
+
+ SetUsePID(kTRUE);
+
+ //activate pileup rejection (for pp)
+ SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ //Do recalculate the vertex
+ SetRemoveDaughtersFromPrim(kTRUE);
+
+ // Cut on the zvtx
+ SetMaxVtxZ(10.);
+
+ // Use the kFirst selection for tracks with candidate D with pt<2
+ SetSelectCandTrackSPDFirst(kTRUE,2.);
+
+ // Use special cuts for D candidates with pt<2
+ SetUseSpecialCuts(kTRUE);
+ SetMaximumPtSpecialCuts(2.);
+
+ PrintAll();
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+}
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
+ //
+ // Final CUTS USED FOR 2010 PbPb analysis of central events
+ // REMEMBER TO EVENTUALLY SWITCH ON
+ // SetUseAOD049(kTRUE);
+ //
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
+
+
+ // CENTRALITY SELECTION
+ SetMinCentrality(0.);
+ SetMaxCentrality(80.);
+ SetUseCentrality(AliRDHFCuts::kCentV0M);
+
+
+
+ // EVENT CUTS
+ SetMinVtxContr(1);
+ // MAX Z-VERTEX CUT
+ SetMaxVtxZ(10.);
+ // PILE UP
+ SetOptPileup(AliRDHFCuts::kNoPileupSelection);
+ // PILE UP REJECTION
+ //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ // TRACKS ON SINGLE TRACKS
+ AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ // esdTrackCuts->SetMinNClustersITS(4);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetPtRange(0.7,1.e10);
+
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+ // SPD k FIRST for Pt<3 GeV/c
+ SetSelectCandTrackSPDFirst(kTRUE, 3);
+
+ // CANDIDATE CUTS
+ const Int_t nptbins =13;
+ const Double_t ptmax = 9999.;
+ const Int_t nvars=11;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=8.;
+ ptbins[9]=12.;
+ ptbins[10]=16.;
+ ptbins[11]=20.;
+ ptbins[12]=24.;
+ ptbins[13]=ptmax;
+
+ SetGlobalIndex(nvars,nptbins);
+ SetPtBins(nptbins+1,ptbins);
+ SetMinPtCandidate(2.);
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.85,0.99,2.},/* pt<0.5*/
+ {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.9,0.99,2.},/* 0.5<pt<1*/
+ {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.995,8.},/* 1<pt<2 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-45000.*1E-8,0.95,0.998,7.},/* 2<pt<3 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
+ {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */
+ {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,8.},/* 12<pt<16 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.995,6.},/* 16<pt<20 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.995,6.},/* 20<pt<24 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.99,2.}};/* pt>24 */
+
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ SetUseSpecialCuts(kTRUE);
+ SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+ cutsMatrixTransposeStand=NULL;
+
+ // PID SETTINGS
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ //pidObj->SetName("pid4D0");
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetPCompatTOF(2.);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+
+
+ SetPidHF(pidObj);
+ SetUsePID(kTRUE);
+ SetUseDefaultPID(kFALSE);
+
+
+ PrintAll();
+
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+
+}
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
+ // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
+
+
+ SetName("D0toKpiCutsStandard");
+ SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
+
+
+ // CENTRALITY SELECTION
+ SetMinCentrality(40.);
+ SetMaxCentrality(80.);
+ SetUseCentrality(AliRDHFCuts::kCentV0M);
+
+ // EVENT CUTS
+ SetMinVtxContr(1);
+ // MAX Z-VERTEX CUT
+ SetMaxVtxZ(10.);
+ // PILE UP
+ SetOptPileup(AliRDHFCuts::kNoPileupSelection);
+ // PILE UP REJECTION
+ //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+ // TRACKS ON SINGLE TRACKS
+ AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+ esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+ esdTrackCuts->SetRequireTPCRefit(kTRUE);
+ esdTrackCuts->SetRequireITSRefit(kTRUE);
+ // esdTrackCuts->SetMinNClustersITS(4);
+ esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+ esdTrackCuts->SetMinDCAToVertexXY(0.);
+ esdTrackCuts->SetEtaRange(-0.8,0.8);
+ esdTrackCuts->SetPtRange(0.5,1.e10);
+
+ esdTrackCuts->SetMaxDCAToVertexXY(1.);
+ esdTrackCuts->SetMaxDCAToVertexZ(1.);
+ esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
+
+
+ AddTrackCuts(esdTrackCuts);
+ delete esdTrackCuts;
+ esdTrackCuts=NULL;
+ // SPD k FIRST for Pt<3 GeV/c
+ SetSelectCandTrackSPDFirst(kTRUE, 3);
+
+ // CANDIDATE CUTS
+ const Int_t nptbins =13;
+ const Double_t ptmax = 9999.;
+ const Int_t nvars=11;
+ Float_t ptbins[nptbins+1];
+ ptbins[0]=0.;
+ ptbins[1]=0.5;
+ ptbins[2]=1.;
+ ptbins[3]=2.;
+ ptbins[4]=3.;
+ ptbins[5]=4.;
+ ptbins[6]=5.;
+ ptbins[7]=6.;
+ ptbins[8]=8.;
+ ptbins[9]=12.;
+ ptbins[10]=16.;
+ ptbins[11]=20.;
+ ptbins[12]=24.;
+ ptbins[13]=ptmax;
+
+ SetGlobalIndex(nvars,nptbins);
+ SetPtBins(nptbins+1,ptbins);
+ SetMinPtCandidate(0.);
+
+ Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-20000.*1E-8,0.7,0.993,2.},/* pt<0.5*/
+ {0.400,400.*1E-4,0.8,0.5,0.5,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.85,0.993,2.},/* 0.5<pt<1*/
+ {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.85,0.995,6.},/* 1<pt<2 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.991,5.},/* 2<pt<3 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.993,5.},/* 3<pt<4 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.995,5.},/* 4<pt<5 */
+ {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
+ {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
+ {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.995,5.},/* 8<pt<12 */
+ {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.83,0.995,4.},/* 12<pt<16 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.82,0.995,4.},/* 16<pt<20 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.81,0.995,4.},/* 20<pt<24 */
+ {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.8,0.995,4.}};/* pt>24 */
+
+
+ //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+ Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+ for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+
+ for (Int_t ibin=0;ibin<nptbins;ibin++){
+ for (Int_t ivar = 0; ivar<nvars; ivar++){
+ cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
+ }
+ }
+
+ SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+ SetUseSpecialCuts(kTRUE);
+ SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
+ for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+ delete [] cutsMatrixTransposeStand;
+ cutsMatrixTransposeStand=NULL;
+
+ // PID SETTINGS
+ AliAODPidHF* pidObj=new AliAODPidHF();
+ //pidObj->SetName("pid4D0");
+ Int_t mode=1;
+ const Int_t nlims=2;
+ Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+ Bool_t compat=kTRUE; //effective only for this mode
+ Bool_t asym=kTRUE;
+ Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+ pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+ pidObj->SetMatch(mode);
+ pidObj->SetPLimit(plims,nlims);
+ pidObj->SetSigma(sigmas);
+ pidObj->SetCompat(compat);
+ pidObj->SetTPC(kTRUE);
+ pidObj->SetTOF(kTRUE);
+ pidObj->SetPCompatTOF(2.);
+ pidObj->SetSigmaForTPCCompat(3.);
+ pidObj->SetSigmaForTOFCompat(3.);
+ pidObj->SetOldPid(kTRUE);
+
+ SetPidHF(pidObj);
+ SetUsePID(kTRUE);
+ SetUseDefaultPID(kFALSE);
+
+ SetLowPt(kTRUE,2.);
+
+ PrintAll();
+
+
+ delete pidObj;
+ pidObj=NULL;
+
+ return;
+
+}
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
+
+ // Default 2010 PbPb cut object
+ SetStandardCutsPbPb2010();
+ AliAODPidHF *pidobj=GetPidHF();
+
+ pidobj->SetOldPid(kFALSE);
+
+ //
+ // Enable all 2011 PbPb run triggers
+ //
+ SetTriggerClass("");
+ ResetMaskAndEnableMBTrigger();
+ EnableCentralTrigger();
+ EnableSemiCentralTrigger();
+
+}
+
+//---------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
+{
+ // ############################################################
+ //
+ // Apply Bayesian PID selection
+ // Implementation of Bayesian PID: Jeremy Wilkinson
+ //
+ // ############################################################
+
+ if (!fUsePID || !d) return 3;
+
+
+ if (fBayesianStrategy == kBayesWeightNoFilter) {
+ //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
+ CalculateBayesianWeights(d);
+ return 3;
+ }
+
+
+
+ Int_t isD0 = 0;
+ Int_t isD0bar = 0;
+ Int_t returnvalue = 0;
+
+ Int_t isPosKaon = 0, isNegKaon = 0, isPosPion = 0, isNegPion = 0;
+
+ //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
+
+
+ Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
+ AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+
+ if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
+ return 0; //reject if charges not opposite
+ }
+ Double_t momentumpositive=0., momentumnegative=0.; //Used in "prob > prior" method
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs
+
+ if (aodtrack[daught]->Charge() == -1) {
+ momentumnegative = aodtrack[daught]->P();
+ }
+ if (aodtrack[daught]->Charge() == +1) {
+ momentumpositive = aodtrack[daught]->P();
+ }
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ checkPIDInfo[daught] = kFALSE;
+ continue;
+ }
+
+ }
+
+ //Loop over daughters ends here
+
+ if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
+ return 0; //Reject if both daughters lack both TPC and TOF info
+ }
+
+
+ CalculateBayesianWeights(d); //Calculates all Bayesian probabilities for both positive and negative tracks
+ // Double_t prob0[AliPID::kSPECIES];
+ // fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
+ ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
+ switch (fBayesianCondition) {
+ ///A: Standard max. probability method (accept most likely species)
+ case kMaxProb:
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
+ isPosKaon = 1; //flag [daught] as a kaon
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kPion]) { //If highest probability lies with pion
+ isPosPion = 1; //flag [daught] as a pion
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
+ isNegKaon = 1; //flag [daught] as a kaon
+ }
+
+ if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kPion]) { //If highest probability lies with kaon
+ isNegPion = 1; //flag [daught] as a pion
+ }
+
+
+ break;
+ ///B: Accept if probability greater than prior
+ case kAbovePrior:
+
+ if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
+ GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) { //Retrieves relevant prior, gets value at momentum
+ isNegKaon = 1;
+ }
+ if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
+ GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) { //Retrieves relevant prior, gets value at momentum
+ isPosKaon = 1;
+ }
+
+ break;
+
+ ///C: Accept if probability greater than user-defined threshold
+ case kThreshold:
+ if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
+ isNegKaon = 1;
+ }
+ if (fWeightsNegative[AliPID::kPion] > fProbThreshold) {
+ isNegPion = 1;
+ }
+
+ if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
+ isPosKaon = 1;
+ }
+
+ if (fWeightsPositive[AliPID::kPion] > fProbThreshold) {
+ isPosPion = 1;
+ }
+
+ break;
+ }
+
+
+ //Momentum-based selection (also applied to filtered weighted method)
+
+ if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
+ if (isNegKaon && isPosKaon) { // If both are kaons, reject
+ isD0 = 1;
+ isD0bar = 1;
+ } else if (isNegKaon && isPosPion) { //If negative kaon present, D0
+ isD0 = 1;
+ } else if (isPosKaon && isNegPion) { //If positive kaon present, D0bar
+ isD0bar = 1;
+ } else { //If neither ID'd as kaon, subject to extra tests
+ isD0 = 1;
+ isD0bar = 1;
+ if (aodtrack[0]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[0], "TOF"))) { //If it's low momentum and not ID'd as a kaon, we assume it must be a pion
+ if (aodtrack[0]->Charge() == -1) {
+ isD0 = 0; //Check charge and reject based on pion hypothesis
+ }
+ if (aodtrack[0]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
+ if (aodtrack[1]->Charge() == -1) {
+ isD0 = 0;
+ }
+ if (aodtrack[1]->Charge() == 1) {
+ isD0bar = 0;
+ }
+ }
+ }
+
+
+
+ if (isD0 && isD0bar) {
+ returnvalue = 3;
+ }
+ if (isD0 && !isD0bar) {
+ returnvalue = 1;
+ }
+ if (!isD0 && isD0bar) {
+ returnvalue = 2;
+ }
+ if (!isD0 && !isD0bar) {
+ returnvalue = 0;
+ }
+ }
+
+ //Simple Bayesian
+ if (fBayesianStrategy == kBayesSimple) {
+
+ if (isPosKaon && isNegKaon) { //If both are ID'd as kaons, accept as possible
+ returnvalue = 3;
+ } else if (isNegKaon && isPosPion) { //If negative kaon, D0
+ returnvalue = 1;
+ } else if (isPosKaon && isNegPion) { //If positive kaon, D0-bar
+ returnvalue = 2;
+ } else if (isPosPion && isNegPion) {
+ returnvalue = 0; //If neither kaon, reject
+ } else {returnvalue = 0;} //default
+
+ }
+
+ return returnvalue;
+
+
+
+}
+
+
+
+//---------------------------------------------------------------------------
+void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
+
+{
+ //Function to compute weights for Bayesian method
+
+ AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
+ if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
+ return; //Reject if charges do not oppose
+ }
+ for (Int_t daught = 0; daught < 2; daught++) {
+ //Loop over prongs
+
+ if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
+ continue;
+ }
+
+
+ // identify kaon, define weights
+ if (aodtrack[daught]->Charge() == +1) {
+ fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
+ }
+
+ if (aodtrack[daught]->Charge() == -1) {
+ fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
+ }
+ }
+}
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff