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