PWGHF/vertexingHF/AliRDHFCutsDstoKKpi.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // Class for cuts on AOD reconstructed Ds->KKpi
  21 //
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TDatabasePDG.h>
  26 #include <Riostream.h>
  27
  28 #include "AliRDHFCutsDstoKKpi.h"
  29 #include "AliAODRecoDecayHF3Prong.h"
  30 #include "AliAODTrack.h"
  31 #include "AliESDtrack.h"
  32
  33 using std::cout;
  34 using std::endl;
  35
  36 ClassImp(AliRDHFCutsDstoKKpi)
  37
  38 //--------------------------------------------------------------------------
  39 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) :
  40 AliRDHFCuts(name),
  41 fPidOption(0),
  42 fMaxPtStrongPid(0.),
  43 fMaxPStrongPidK(0.),
  44 fMaxPStrongPidpi(0.)
  45 {
  46   //
  47   // Default Constructor
  48   //
  49   Int_t nvars=20;
  50   SetNVars(nvars);
  51   TString varNames[20]={"inv. mass [GeV]",
  52                         "pTK [GeV/c]",
  53                         "pTPi [GeV/c]",
  54                         "d0K [cm]",
  55                         "d0Pi [cm]",
  56                         "dist12 [cm]",
  57                         "sigmavert [cm]",
  58                         "decLen [cm]",
  59                         "ptMax [GeV/c]",
  60                         "cosThetaPoint",
  61                         "Sum d0^2 (cm^2)",
  62                         "dca [cm]",
  63                         "inv. mass (Mphi-MKK) [GeV]",
  64                         "inv. mass (MKo*-MKpi) [GeV]",
  65                         "Abs(CosineKpiPhiRFrame)^3",
  66                         "CosPiDsLabFrame",
  67                         "decLenXY [cm]"
  68                         "NormdecLen",
  69                         "NormdecLenXY [cm]",
  70                         "cosThetaPointXY"};
  71
  72   Bool_t isUpperCut[20]={kTRUE,
  73                          kFALSE,
  74                          kFALSE,
  75                          kFALSE,
  76                          kFALSE,
  77                          kFALSE,
  78                          kTRUE,
  79                          kFALSE,
  80                          kFALSE,
  81                          kFALSE,
  82                          kFALSE,
  83                          kTRUE,
  84                          kTRUE,
  85                          kTRUE,
  86                          kFALSE,
  87                          kTRUE,
  88                          kFALSE,
  89                          kFALSE,
  90                          kFALSE,
  91                          kFALSE};
  92   SetVarNames(20,varNames,isUpperCut);
  93   Bool_t forOpt[20]={kFALSE,
  94                     kFALSE,
  95                     kFALSE,
  96                     kFALSE,
  97                     kFALSE,
  98                     kFALSE,
  99                     kTRUE,
 100                     kTRUE,
 101                     kTRUE,
 102                     kTRUE,
 103                     kTRUE,
 104                     kFALSE,
 105                     kTRUE,
 106                     kTRUE,
 107                     kFALSE,
 108                     kFALSE,
 109                     kTRUE,
 110                     kTRUE,
 111                     kTRUE,
 112                     kTRUE};
 113
 114   SetVarsForOpt(11,forOpt);
 115   Float_t limits[2]={0,999999999.};
 116   SetPtBins(2,limits);
 117   if(fPidHF)delete fPidHF;
 118   fPidHF=new AliAODPidHF();
 119   Double_t plim[2]={0.6,0.8};
 120   Double_t nsigma[5]={2.,1.,2.,3.,0.};
 121
 122   fPidHF->SetPLimit(plim);
 123   fPidHF->SetAsym(kTRUE);
 124   fPidHF->SetSigma(nsigma);
 125   fPidHF->SetMatch(1);
 126   fPidHF->SetTPC(1);
 127   fPidHF->SetTOF(1);
 128   fPidHF->SetITS(0);
 129   fPidHF->SetTRD(0);
 130   fPidHF->SetCompat(kTRUE);
 131
 132 }
 133 //--------------------------------------------------------------------------
 134 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
 135   AliRDHFCuts(source),
 136   fPidOption(source.fPidOption),
 137   fMaxPtStrongPid(source.fMaxPtStrongPid),
 138   fMaxPStrongPidK(source.fMaxPStrongPidK),
 139   fMaxPStrongPidpi(source.fMaxPStrongPidpi)
 140 {
 141   //
 142   // Copy constructor
 143   //
 144
 145 }
 146 //--------------------------------------------------------------------------
 147 AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
 148 {
 149   //
 150   // assignment operator
 151   //
 152   if(&source == this) return *this;
 153
 154   AliRDHFCuts::operator=(source);
 155
 156   fPidOption=source.fPidOption;
 157   fMaxPtStrongPid=source.fMaxPtStrongPid;
 158   fMaxPStrongPidK=source.fMaxPStrongPidK;
 159   fMaxPStrongPidpi=source.fMaxPStrongPidpi;
 160
 161   return *this;
 162 }
 163
 164
 165 //---------------------------------------------------------------------------
 166 void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
 167   //
 168   // Fills in vars the values of the variables
 169   //
 170
 171   if(nvars!=fnVarsForOpt) {
 172     printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
 173     return;
 174   }
 175
 176   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 177
 178   //recalculate vertex w/o daughters
 179   Bool_t cleanvtx=kFALSE;
 180   AliAODVertex *origownvtx=0x0;
 181   if(fRemoveDaughtersFromPrimary) {
 182     if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
 183     cleanvtx=kTRUE;
 184     if(!RecalcOwnPrimaryVtx(dd,aod)) {
 185       CleanOwnPrimaryVtx(dd,aod,origownvtx);
 186       cleanvtx=kFALSE;
 187     }
 188   }
 189
 190   Int_t iter=-1;
 191   if(fVarsForOpt[0]){
 192     iter++;
 193     if(TMath::Abs(pdgdaughters[0])==321){
 194       vars[iter]=dd->InvMassDsKKpi();
 195     }else{
 196       vars[iter]=dd->InvMassDspiKK();
 197     }
 198   }
 199   if(fVarsForOpt[1]){
 200     iter++;
 201     Float_t minPtDau=99999.;
 202     for(Int_t iprong=0;iprong<3;iprong++){
 203       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 204          dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
 205     }
 206     vars[iter]=minPtDau;
 207   }
 208   if(fVarsForOpt[2]){
 209     iter++;
 210     for(Int_t iprong=0;iprong<3;iprong++){
 211       if(TMath::Abs(pdgdaughters[iprong])==211) {
 212         vars[iter]=dd->PtProng(iprong);
 213       }
 214     }
 215   }
 216   if(fVarsForOpt[3]){
 217     iter++;
 218     Float_t minImpParDau=99999.;
 219     for(Int_t iprong=0;iprong<3;iprong++){
 220       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 221          dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
 222     }
 223     vars[iter]=minImpParDau;
 224   }
 225   if(fVarsForOpt[4]){
 226     iter++;
 227     for(Int_t iprong=0;iprong<3;iprong++){
 228       if(TMath::Abs(pdgdaughters[iprong])==211) {
 229         vars[iter]=dd->Getd0Prong(iprong);
 230       }
 231     }
 232   }
 233   if(fVarsForOpt[5]){
 234     iter++;
 235     Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
 236     vars[iter]=minDistPair;
 237   }
 238   if(fVarsForOpt[6]){
 239     iter++;
 240     vars[iter]=dd->GetSigmaVert(aod);
 241   }
 242   if(fVarsForOpt[7]){
 243     iter++;
 244     vars[iter] = dd->DecayLength();
 245   }
 246   if(fVarsForOpt[8]){
 247     iter++;
 248     Float_t ptmax=0;
 249     for(Int_t i=0;i<3;i++){
 250       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
 251     }
 252     vars[iter]=ptmax;
 253   }
 254   if(fVarsForOpt[9]){
 255     iter++;
 256     vars[iter]=dd->CosPointingAngle();
 257   }
 258   if(fVarsForOpt[10]){
 259     iter++;
 260     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
 261   }
 262   if(fVarsForOpt[11]){
 263     iter++;
 264     Float_t maxDCA=0.;
 265     for(Int_t i=0;i<3;i++){
 266       if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
 267     }
 268     vars[iter]=maxDCA;
 269   }
 270   if(fVarsForOpt[12]){
 271     iter++;
 272     Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 273     if(TMath::Abs(pdgdaughters[0])==321){
 274
 275       Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
 276        vars[iter]=TMath::Abs(phimass01-mPDGPhi);
 277        // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
 278     }else{
 279       Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
 280        vars[iter]=TMath::Abs(phimass12-mPDGPhi);
 281        // vars[iter]=dd->InvMass2Prongs(1,2,321,321);
 282     }
 283   }
 284   if(fVarsForOpt[13]){
 285     iter++;
 286     Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 287     if(TMath::Abs(pdgdaughters[0])==321){
 288
 289       Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
 290       vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
 291       //              vars[iter]=dd->InvMass2Prongs(1,2,321,211);
 292     }else{
 293       Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
 294       vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
 295       //        vars[iter]=dd->InvMass2Prongs(0,1,211,321);
 296     }
 297   }
 298   if(fVarsForOpt[14]){
 299     iter++;
 300     if(TMath::Abs(pdgdaughters[0])==321){
 301       vars[iter]=dd->CosPiKPhiRFrameKKpi();
 302     }else{
 303       vars[iter]=dd->CosPiKPhiRFramepiKK();
 304     }
 305   }
 306   if(fVarsForOpt[15]){
 307     iter++;
 308     if(TMath::Abs(pdgdaughters[0])==321){
 309       vars[iter]=dd->CosPiDsLabFrameKKpi();
 310     }else{
 311       vars[iter]=dd->CosPiDsLabFramepiKK();
 312     }
 313   }
 314
 315   if(fVarsForOpt[16]){
 316     iter++;
 317     vars[iter]=dd->DecayLengthXY();
 318   }
 319
 320   if(fVarsForOpt[17]){
 321     iter++;
 322     vars[iter]=dd->NormalizedDecayLength();
 323   }
 324
 325   if(fVarsForOpt[18]){
 326     iter++;
 327     vars[iter]=dd->NormalizedDecayLengthXY();
 328   }
 329
 330   if(fVarsForOpt[19]){
 331     iter++;
 332     vars[iter]=dd->CosPointingAngleXY();
 333   }
 334
 335   if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
 336   return;
 337 }
 338 //---------------------------------------------------------------------------
 339 Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 340 {
 341   //
 342   // Checking if Ds is in fiducial acceptance region
 343   //
 344
 345   if(pt > 5.) {
 346     // applying cut for pt > 5 GeV
 347     AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt));
 348     if (TMath::Abs(y) > 0.8) return kFALSE;
 349
 350   } else {
 351     // appliying smooth cut for pt < 5 GeV
 352     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 353     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 354     AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 355     if (y < minFiducialY || y > maxFiducialY) return kFALSE;
 356   }
 357
 358   return kTRUE;
 359 }
 360
 361 //---------------------------------------------------------------------------
 362 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
 363   // PID selection
 364   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both
 365   Int_t retCode=3;
 366   Bool_t okKKpi=kTRUE;
 367   Bool_t okpiKK=kTRUE;
 368   if(!fUsePID || !rd) return retCode;
 369   if(!fPidHF){
 370     AliWarning("AliAODPidHF not created!");
 371     return retCode;
 372   }
 373
 374   Double_t origCompatTOF=fPidHF->GetPCompatTOF();
 375   Double_t origThreshTPC=fPidHF->GetPtThresholdTPC();
 376   if(fPidOption==kStrong){
 377     fPidHF->SetPCompatTOF(999999.);
 378     fPidHF->SetPtThresholdTPC(999999.);
 379   }
 380
 381   Int_t nKaons=0;
 382   Int_t nNotKaons=0;
 383   Int_t sign= rd->GetCharge();
 384   for(Int_t iDaught=0; iDaught<3; iDaught++){
 385     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
 386
 387     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 388     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 389     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
 390
 391     if(isProton>0 &&  isKaon<0  && isPion<0){
 392       fPidHF->SetPCompatTOF(origCompatTOF);
 393       fPidHF->SetPtThresholdTPC(origThreshTPC);
 394       return 0;
 395     }
 396     if(sign!=track->Charge()){// must be kaon
 397       if(isKaon<0){
 398         fPidHF->SetPCompatTOF(origCompatTOF);
 399         fPidHF->SetPtThresholdTPC(origThreshTPC);
 400         return 0;
 401       }
 402       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid && isKaon<=0){
 403         fPidHF->SetPCompatTOF(origCompatTOF);
 404         fPidHF->SetPtThresholdTPC(origThreshTPC);
 405         return 0;
 406       }
 407       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 408         if(isKaon<=0 && track->P()<fMaxPStrongPidK) return 0;
 409       }
 410     }
 411
 412     if(isKaon>0 && isPion<0) nKaons++;
 413     if(isKaon<0) nNotKaons++;
 414     if(iDaught==0){
 415       if(isKaon<0) okKKpi=kFALSE;
 416       if(isPion<0) okpiKK=kFALSE;
 417       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
 418         if(isKaon<=0) okKKpi=kFALSE;
 419         if(isPion<=0) okpiKK=kFALSE;
 420       }
 421       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 422         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okKKpi=kFALSE;
 423         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okpiKK=kFALSE;
 424       }
 425     }
 426     else if(iDaught==2){
 427       if(isKaon<0) okpiKK=kFALSE;
 428       if(isPion<0) okKKpi=kFALSE;
 429        if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
 430          if(isKaon<=0) okpiKK=kFALSE;
 431          if(isPion<=0) okKKpi=kFALSE;
 432       }
 433       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 434         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okpiKK=kFALSE;
 435         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okKKpi=kFALSE;
 436       }
 437     }
 438   }
 439
 440   fPidHF->SetPCompatTOF(origCompatTOF);
 441   fPidHF->SetPtThresholdTPC(origThreshTPC);
 442
 443   if(nKaons>2)return 0;
 444   if(nNotKaons>1) return 0;
 445
 446   if(!okKKpi) retCode-=1;
 447   if(!okpiKK) retCode-=2;
 448
 449   return retCode;
 450 }
 451
 452 //---------------------------------------------------------------------------
 453 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
 454   //
 455   // Apply selection
 456   //
 457
 458   if(!fCutsRD){
 459     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 460     return 0;
 461   }
 462   //PrintAll();
 463   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 464
 465   if(!d){
 466     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
 467     return 0;
 468   }
 469
 470   if(fKeepSignalMC) if(IsSignalMC(d,aod,431)) return 3;
 471
 472   Double_t ptD=d->Pt();
 473   if(ptD<fMinPtCand) return 0;
 474   if(ptD>fMaxPtCand) return 0;
 475
 476   if(d->HasBadDaughters()) return 0;
 477
 478
 479   // selection on daughter tracks
 480   if(selectionLevel==AliRDHFCuts::kAll ||
 481      selectionLevel==AliRDHFCuts::kTracks) {
 482     if(!AreDaughtersSelected(d)) return 0;
 483   }
 484
 485
 486
 487
 488   // selection on candidate
 489   if(selectionLevel==AliRDHFCuts::kAll ||
 490      selectionLevel==AliRDHFCuts::kCandidate) {
 491     //recalculate vertex w/o daughters
 492     AliAODVertex *origownvtx=0x0;
 493     if(fRemoveDaughtersFromPrimary) {
 494       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 495       if(!RecalcOwnPrimaryVtx(d,aod)) {
 496         CleanOwnPrimaryVtx(d,aod,origownvtx);
 497         return 0;
 498       }
 499     }
 500
 501     Int_t okDsKKpi=1;
 502     Int_t okDspiKK=1;
 503     Int_t okMassPhiKKpi=0;
 504     Int_t okMassPhipiKK=0;
 505     Int_t okMassK0starKKpi=0;
 506     Int_t okMassK0starpiKK=0;
 507     Int_t okDsPhiKKpi=0;
 508     Int_t okDsPhipiKK=0;
 509     Int_t okDsK0starKKpi=0;
 510     Int_t okDsK0starpiKK=0;
 511
 512     Double_t pt=d->Pt();
 513     Int_t ptbin=PtBin(pt);
 514     if (ptbin==-1) {
 515       CleanOwnPrimaryVtx(d,aod,origownvtx);
 516       return 0;
 517     }
 518
 519     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
 520     Double_t mDsKKpi=d->InvMassDsKKpi();
 521     Double_t mDspiKK=d->InvMassDspiKK();
 522     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
 523     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
 524     if(!okDsKKpi && !okDspiKK){
 525       CleanOwnPrimaryVtx(d,aod,origownvtx);
 526       return 0;
 527     }
 528
 529
 530
 531     // cuts on resonant decays (via Phi or K0*)
 532     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 533     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 534     if(okDsKKpi){
 535       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
 536       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
 537       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhiKKpi=1;
 538       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starKKpi = 1;
 539       if(!okMassPhiKKpi && !okMassK0starKKpi) okDsKKpi=0;
 540       if(okMassPhiKKpi) okDsPhiKKpi=1;
 541       if(okMassK0starKKpi) okDsK0starKKpi=1;
 542     }
 543     if(okDspiKK){
 544       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
 545       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
 546       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starpiKK = 1;
 547       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhipiKK=1;
 548       if(!okMassPhipiKK && !okMassK0starpiKK) okDspiKK=0;
 549       if(okMassPhipiKK) okDsPhipiKK=1;
 550       if(okMassK0starpiKK) okDsK0starpiKK=1;
 551     }
 552     if(!okDsKKpi && !okDspiKK){
 553       CleanOwnPrimaryVtx(d,aod,origownvtx);
 554       return 0;
 555     }
 556
 557     // Cuts on track pairs
 558     for(Int_t i=0;i<3;i++){
 559       if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]){
 560         CleanOwnPrimaryVtx(d,aod,origownvtx);
 561         return 0;
 562       }
 563     }
 564     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] ||
 565        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]){
 566       CleanOwnPrimaryVtx(d,aod,origownvtx);
 567       return 0;
 568     }
 569
 570
 571
 572     //single track
 573     if(TMath::Abs(d->Pt2Prong(1)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 574        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]){
 575       CleanOwnPrimaryVtx(d,aod,origownvtx);
 576       return 0;
 577     }
 578
 579     if(okDsKKpi){
 580       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 581          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
 582       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 583          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
 584     }
 585     if(okDspiKK){
 586       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 587          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
 588       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 589          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
 590     }
 591     if(!okDsKKpi && !okDspiKK){
 592       CleanOwnPrimaryVtx(d,aod,origownvtx);
 593       return 0;
 594     }
 595
 596     // Cuts on candidate triplet
 597
 598
 599     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]){
 600       CleanOwnPrimaryVtx(d,aod,origownvtx);
 601       return 0;
 602     }
 603
 604     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 605        d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 606        d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {
 607       CleanOwnPrimaryVtx(d,aod,origownvtx);
 608       return 0;
 609     }
 610
 611     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]){
 612       CleanOwnPrimaryVtx(d,aod,origownvtx);
 613       return 0;
 614     }
 615
 616
 617     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 618     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]){
 619       CleanOwnPrimaryVtx(d,aod,origownvtx);
 620       return 0;
 621     }
 622
 623
 624     //sec vert
 625     Double_t sigmavert=d->GetSigmaVert(aod);
 626     if(sigmavert>fCutsRD[GetGlobalIndex(6,ptbin)]){
 627       CleanOwnPrimaryVtx(d,aod,origownvtx);
 628       return 0;
 629     }
 630
 631     // decay length XY
 632     if(d->DecayLengthXY()<fCutsRD[GetGlobalIndex(16,ptbin)]){
 633       CleanOwnPrimaryVtx(d,aod,origownvtx);
 634       return 0;
 635     }
 636
 637     //norm decay length
 638     if(d->NormalizedDecayLength()<fCutsRD[GetGlobalIndex(17,ptbin)]){
 639       CleanOwnPrimaryVtx(d,aod,origownvtx);
 640       return 0;
 641     }
 642
 643     //norm decay length XY
 644     if(d->NormalizedDecayLengthXY()<fCutsRD[GetGlobalIndex(18,ptbin)]){
 645       CleanOwnPrimaryVtx(d,aod,origownvtx);
 646       return 0;
 647     }
 648
 649     //cos pointing XY
 650     if(d->CosPointingAngleXY()<fCutsRD[GetGlobalIndex(19,ptbin)]){
 651       CleanOwnPrimaryVtx(d,aod,origownvtx);
 652       return 0;
 653     }
 654
 655
 656     if(okDsKKpi){
 657       Double_t cosPiKPhiRFKKpi=d->CosPiKPhiRFrameKKpi();
 658       Double_t kincutPiKPhiKKpi=TMath::Abs(cosPiKPhiRFKKpi*cosPiKPhiRFKKpi*cosPiKPhiRFKKpi);
 659       if(kincutPiKPhiKKpi<fCutsRD[GetGlobalIndex(14,ptbin)]) okDsKKpi=0;
 660     }
 661     if(okDspiKK){
 662       Double_t cosPiKPhiRFpiKK=d->CosPiKPhiRFramepiKK();
 663       Double_t kincutPiKPhipiKK=TMath::Abs(cosPiKPhiRFpiKK*cosPiKPhiRFpiKK*cosPiKPhiRFpiKK);
 664       if(kincutPiKPhipiKK<fCutsRD[GetGlobalIndex(14,ptbin)]) okDspiKK=0;
 665     }
 666     if(!okDsKKpi && !okDspiKK){
 667       CleanOwnPrimaryVtx(d,aod,origownvtx);
 668       return 0;
 669     }
 670
 671
 672
 673     if(okDsKKpi){
 674       Double_t cosPiDsLabFrameKKpi=d->CosPiDsLabFrameKKpi();
 675       if(cosPiDsLabFrameKKpi>fCutsRD[GetGlobalIndex(15,ptbin)]) okDsKKpi=0;
 676     }
 677     if(okDspiKK){
 678       Double_t cosPiDsLabFramepiKK=d->CosPiDsLabFramepiKK();
 679       if(cosPiDsLabFramepiKK>fCutsRD[GetGlobalIndex(15,ptbin)]) okDspiKK=0;
 680     }
 681     if(!okDsKKpi && !okDspiKK){
 682       CleanOwnPrimaryVtx(d,aod,origownvtx);
 683       return 0;
 684     }
 685
 686      // unset recalculated primary vertex when not needed any more
 687     CleanOwnPrimaryVtx(d,aod,origownvtx);
 688
 689
 690
 691     if(!okDsKKpi){
 692       okDsPhiKKpi=0;
 693       okDsK0starKKpi=0;
 694     }
 695     if(!okDspiKK){
 696       okDsPhipiKK=0;
 697       okDsK0starpiKK=0;
 698     }
 699
 700     // PID selection
 701     Int_t returnvaluePID=3;
 702     if(selectionLevel==AliRDHFCuts::kAll ||
 703        selectionLevel==AliRDHFCuts::kCandidate ||
 704        selectionLevel==AliRDHFCuts::kPID) {
 705       returnvaluePID = IsSelectedPID(d);
 706       fIsSelectedPID=returnvaluePID;
 707     }
 708     if(returnvaluePID==0)return 0;
 709
 710     Bool_t okPidDsKKpi=returnvaluePID&1;
 711     Bool_t okPidDspiKK=returnvaluePID&2;
 712     if(!okPidDsKKpi){
 713       okDsPhiKKpi=0;
 714       okDsK0starKKpi=0;
 715     }
 716     if(!okPidDspiKK){
 717       okDsPhipiKK=0;
 718       okDsK0starpiKK=0;
 719     }
 720
 721     if((okPidDsKKpi && okDsKKpi)||(okPidDspiKK && okDspiKK)){
 722       Int_t returnvalue=0;
 723       if(okDsKKpi) returnvalue+=1;
 724       if(okDspiKK) returnvalue+=2;
 725       if(okDsPhiKKpi) returnvalue+=4;
 726       if(okDsPhipiKK) returnvalue+=8;
 727       if(okDsK0starKKpi) returnvalue+=16;
 728       if(okDsK0starpiKK) returnvalue+=32;
 729       return returnvalue;
 730     }else{
 731       return 0;
 732     }
 733   }
 734   return 15;
 735
 736 }
 737
 738 //--------------------------------------------------------------------------
 739
 740 UInt_t AliRDHFCutsDstoKKpi::GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const{
 741
 742   UInt_t bitmap=0;
 743
 744   Double_t sigmaTPCPionHyp=-999.;
 745   Double_t sigmaTPCKaonHyp=-999.;
 746   Double_t sigmaTPCProtonHyp=-999.;
 747   Double_t sigmaTOFPionHyp=-999.;
 748   Double_t sigmaTOFKaonHyp=-999.;
 749   Double_t sigmaTOFProtonHyp=-999.;
 750
 751   Int_t oksigmaTPCPionHyp=fPidHF->GetnSigmaTPC(track,2,sigmaTPCPionHyp);
 752   Int_t oksigmaTPCKaonHyp=fPidHF->GetnSigmaTPC(track,3,sigmaTPCKaonHyp);
 753   Int_t oksigmaTPCProtonHyp=fPidHF->GetnSigmaTPC(track,4,sigmaTPCProtonHyp);
 754   Int_t oksigmaTOFPionHyp=fPidHF->GetnSigmaTOF(track,2,sigmaTOFPionHyp);
 755   Int_t oksigmaTOFKaonHyp=fPidHF->GetnSigmaTOF(track,3,sigmaTOFKaonHyp);
 756   Int_t oksigmaTOFProtonHyp=fPidHF->GetnSigmaTOF(track,4,sigmaTOFProtonHyp);
 757
 758   sigmaTPCPionHyp=TMath::Abs(sigmaTPCPionHyp);
 759   sigmaTPCKaonHyp=TMath::Abs(sigmaTPCKaonHyp);
 760   sigmaTPCProtonHyp=TMath::Abs(sigmaTPCProtonHyp);
 761   sigmaTOFPionHyp=TMath::Abs(sigmaTOFPionHyp);
 762   sigmaTOFKaonHyp=TMath::Abs(sigmaTOFKaonHyp);
 763   sigmaTOFProtonHyp=TMath::Abs(sigmaTOFProtonHyp);
 764
 765  if (oksigmaTPCPionHyp && sigmaTPCPionHyp>0.){
 766     if (sigmaTPCPionHyp<1.) bitmap+=1<<kTPCPionLess1;
 767     else{
 768       if (sigmaTPCPionHyp<2.) bitmap+=1<<kTPCPionMore1Less2;
 769       else {
 770         if (sigmaTPCPionHyp<3.) bitmap+=1<<kTPCPionMore2Less3;
 771         else bitmap+=1<<kTPCPionMore3;
 772       }
 773     }
 774   }
 775
 776   if (oksigmaTPCKaonHyp && sigmaTPCKaonHyp>0.){
 777     if (sigmaTPCKaonHyp<1.) bitmap+=1<<kTPCKaonLess1;
 778     else{
 779       if (sigmaTPCKaonHyp<2.) bitmap+=1<<kTPCKaonMore1Less2;
 780       else {
 781         if (sigmaTPCKaonHyp<3.) bitmap+=1<<kTPCKaonMore2Less3;
 782         else bitmap+=1<<kTPCKaonMore3;
 783       }
 784     }
 785   }
 786
 787   if (oksigmaTPCProtonHyp && sigmaTPCProtonHyp>0.){
 788     if (sigmaTPCProtonHyp<1.) bitmap+=1<<kTPCProtonLess1;
 789     else{
 790       if (sigmaTPCProtonHyp<2.) bitmap+=1<<kTPCProtonMore1Less2;
 791       else {
 792         if (sigmaTPCProtonHyp<3.) bitmap+=1<<kTPCProtonMore2Less3;
 793         else bitmap+=1<<kTPCProtonMore3;
 794       }
 795     }
 796   }
 797
 798   if (oksigmaTOFPionHyp && sigmaTOFPionHyp>0.){
 799     if (sigmaTOFPionHyp<1.) bitmap+=1<<kTOFPionLess1;
 800     else{
 801       if (sigmaTOFPionHyp<2.) bitmap+=1<<kTOFPionMore1Less2;
 802       else {
 803         if (sigmaTOFPionHyp<3.) bitmap+=1<<kTOFPionMore2Less3;
 804         else bitmap+=1<<kTOFPionMore3;
 805       }
 806     }
 807   }
 808
 809   if (oksigmaTOFKaonHyp && sigmaTOFKaonHyp>0.){
 810     if (sigmaTOFKaonHyp<1.) bitmap+=1<<kTOFKaonLess1;
 811     else{
 812       if (sigmaTOFKaonHyp<2.) bitmap+=1<<kTOFKaonMore1Less2;
 813       else {
 814         if (sigmaTOFKaonHyp<3.) bitmap+=1<<kTOFKaonMore2Less3;
 815         else bitmap+=1<<kTOFKaonMore3;
 816       }
 817     }
 818   }
 819
 820   if (oksigmaTOFProtonHyp && sigmaTOFProtonHyp>0.){
 821     if (sigmaTOFProtonHyp<1.) bitmap+=1<<kTOFProtonLess1;
 822     else{
 823       if (sigmaTOFProtonHyp<2.) bitmap+=1<<kTOFProtonMore1Less2;
 824       else {
 825         if (sigmaTOFProtonHyp<3.) bitmap+=1<<kTOFProtonMore2Less3;
 826         else bitmap+=1<<kTOFProtonMore3;
 827       }
 828     }
 829   }
 830
 831
 832
 833   return bitmap;
 834
 835 }
 836
 837
 838 //---------------------------------------------------------------------------
 839
 840 void AliRDHFCutsDstoKKpi::SetStandardCutsPP2010() {
 841   //
 842   //STANDARD CUTS USED FOR 2010 pp analysis
 843   //
 844
 845   SetName("DstoKKpiCutsStandard");
 846   SetTitle("Standard Cuts for D+s analysis");
 847
 848   // PILE UP REJECTION
 849   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 850
 851   // EVENT CUTS
 852   SetMinVtxContr(1);
 853
 854   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 855   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 856   //default
 857   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 858   esdTrackCuts->SetRequireITSRefit(kTRUE);
 859   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
 860   esdTrackCuts->SetMinNClustersTPC(70);
 861   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 862                                          AliESDtrackCuts::kAny);
 863   // default is kBoth, otherwise kAny
 864   esdTrackCuts->SetMinDCAToVertexXY(0.);
 865   esdTrackCuts->SetPtRange(0.3,1.e10);
 866
 867   AddTrackCuts(esdTrackCuts);
 868
 869
 870
 871   const Int_t nptbins=4;
 872   Float_t ptbins[nptbins+1];
 873   ptbins[0]=2.;
 874   ptbins[1]=4.;
 875   ptbins[2]=6.;
 876   ptbins[3]=8.;
 877   ptbins[4]=12.;
 878
 879   const Int_t nvars=20;
 880
 881   Float_t** anacutsval;
 882   anacutsval=new Float_t*[nvars];
 883
 884   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
 885   for(Int_t ipt=0;ipt<nptbins;ipt++){
 886
 887     anacutsval[0][ipt]=0.35;
 888     anacutsval[1][ipt]=0.3;
 889     anacutsval[2][ipt]=0.3;
 890     anacutsval[3][ipt]=0.;
 891     anacutsval[4][ipt]=0.;
 892     anacutsval[5][ipt]=0.005;
 893     anacutsval[8][ipt]=0.;
 894     anacutsval[10][ipt]=0.;
 895     anacutsval[11][ipt]=1000.0;
 896     anacutsval[13][ipt]=0.1;
 897     anacutsval[16][ipt]=0.;
 898     anacutsval[17][ipt]=0.;
 899     anacutsval[18][ipt]=0.;
 900     anacutsval[19][ipt]=-1.;
 901
 902
 903   }
 904
 905   //sigmavertex
 906
 907   anacutsval[6][0]=0.020;
 908   anacutsval[6][1]=0.030;
 909   anacutsval[6][2]=0.030;
 910   anacutsval[6][3]=0.060;
 911
 912   //Decay length
 913
 914   anacutsval[7][0]=0.035;
 915   anacutsval[7][1]=0.035;
 916   anacutsval[7][2]=0.040;
 917   anacutsval[7][3]=0.040;
 918
 919   //cosThetaPoint
 920
 921   anacutsval[9][0]=0.94;
 922   anacutsval[9][1]=0.94;
 923   anacutsval[9][2]=0.94;
 924   anacutsval[9][3]=0.94;
 925
 926   //phi DeltaMass
 927
 928   anacutsval[12][0]=0.0080;
 929   anacutsval[12][1]=0.0050;
 930   anacutsval[12][2]=0.0045;
 931   anacutsval[12][3]=0.0090;
 932
 933   //Kin1
 934
 935   anacutsval[14][0]=0.10;
 936   anacutsval[14][1]=0.05;
 937   anacutsval[14][2]=0.0;
 938   anacutsval[14][3]=0.05;
 939
 940   //Kin2
 941
 942   anacutsval[15][0]=0.95;
 943   anacutsval[15][1]=0.95;
 944   anacutsval[15][2]=1.;
 945   anacutsval[15][3]=0.95;
 946
 947
 948   SetUsePID(kTRUE);
 949   SetPidOption(1);
 950   SetMaxPtStrongPid(9999.);
 951   SetGlobalIndex(nvars,nptbins);
 952   SetPtBins(nptbins+1,ptbins);
 953   SetCuts(nvars,nptbins,anacutsval);
 954   SetRemoveDaughtersFromPrim(kTRUE);
 955
 956   PrintAll();
 957
 958   for(Int_t iic=0;iic<nvars;iic++){delete [] anacutsval[iic];}
 959   delete [] anacutsval;
 960   anacutsval=NULL;
 961
 962   return;
 963 }
 964