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(fMaxRapidityCand>-998.){
 346     if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
 347     else return kTRUE;
 348   }
 349
 350   if(pt > 5.) {
 351     // applying cut for pt > 5 GeV
 352     AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt));
 353     if (TMath::Abs(y) > 0.8) return kFALSE;
 354
 355   } else {
 356     // appliying smooth cut for pt < 5 GeV
 357     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 358     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 359     AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 360     if (y < minFiducialY || y > maxFiducialY) return kFALSE;
 361   }
 362
 363   return kTRUE;
 364 }
 365
 366 //---------------------------------------------------------------------------
 367 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
 368   // PID selection
 369   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both
 370   Int_t retCode=3;
 371   Bool_t okKKpi=kTRUE;
 372   Bool_t okpiKK=kTRUE;
 373   if(!fUsePID || !rd) return retCode;
 374   if(!fPidHF){
 375     AliWarning("AliAODPidHF not created!");
 376     return retCode;
 377   }
 378
 379   Double_t origCompatTOF=fPidHF->GetPCompatTOF();
 380   Double_t origThreshTPC=fPidHF->GetPtThresholdTPC();
 381   if(fPidOption==kStrong){
 382     fPidHF->SetPCompatTOF(999999.);
 383     fPidHF->SetPtThresholdTPC(999999.);
 384   }
 385
 386   Int_t nKaons=0;
 387   Int_t nNotKaons=0;
 388   Int_t sign= rd->GetCharge();
 389   for(Int_t iDaught=0; iDaught<3; iDaught++){
 390     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
 391
 392     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 393     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 394     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
 395
 396     if(isProton>0 &&  isKaon<0  && isPion<0){
 397       fPidHF->SetPCompatTOF(origCompatTOF);
 398       fPidHF->SetPtThresholdTPC(origThreshTPC);
 399       return 0;
 400     }
 401     if(sign!=track->Charge()){// must be kaon
 402       if(isKaon<0){
 403         fPidHF->SetPCompatTOF(origCompatTOF);
 404         fPidHF->SetPtThresholdTPC(origThreshTPC);
 405         return 0;
 406       }
 407       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid && isKaon<=0){
 408         fPidHF->SetPCompatTOF(origCompatTOF);
 409         fPidHF->SetPtThresholdTPC(origThreshTPC);
 410         return 0;
 411       }
 412       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 413         if(isKaon<=0 && track->P()<fMaxPStrongPidK) return 0;
 414       }
 415     }
 416
 417     if(isKaon>0 && isPion<0) nKaons++;
 418     if(isKaon<0) nNotKaons++;
 419     if(iDaught==0){
 420       if(isKaon<0) okKKpi=kFALSE;
 421       if(isPion<0) okpiKK=kFALSE;
 422       if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
 423         if(isKaon<=0) okKKpi=kFALSE;
 424         if(isPion<=0) okpiKK=kFALSE;
 425       }
 426       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 427         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okKKpi=kFALSE;
 428         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okpiKK=kFALSE;
 429       }
 430     }
 431     else if(iDaught==2){
 432       if(isKaon<0) okpiKK=kFALSE;
 433       if(isPion<0) okKKpi=kFALSE;
 434        if(fPidOption==kStrong && rd->Pt()<fMaxPtStrongPid){
 435          if(isKaon<=0) okpiKK=kFALSE;
 436          if(isPion<=0) okKKpi=kFALSE;
 437       }
 438       if(fPidOption==kStrongPDep && rd->Pt()<fMaxPtStrongPid){
 439         if(isKaon<=0 && track->P()<fMaxPStrongPidK) okpiKK=kFALSE;
 440         if(isPion<=0 && track->P()<fMaxPStrongPidpi) okKKpi=kFALSE;
 441       }
 442     }
 443   }
 444
 445   fPidHF->SetPCompatTOF(origCompatTOF);
 446   fPidHF->SetPtThresholdTPC(origThreshTPC);
 447
 448   if(nKaons>2)return 0;
 449   if(nNotKaons>1) return 0;
 450
 451   if(!okKKpi) retCode-=1;
 452   if(!okpiKK) retCode-=2;
 453
 454   return retCode;
 455 }
 456
 457 //---------------------------------------------------------------------------
 458 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
 459   //
 460   // Apply selection
 461   //
 462
 463   if(!fCutsRD){
 464     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 465     return 0;
 466   }
 467   //PrintAll();
 468   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 469
 470   if(!d){
 471     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
 472     return 0;
 473   }
 474
 475   if(fKeepSignalMC) if(IsSignalMC(d,aod,431)) return 3;
 476
 477   Double_t ptD=d->Pt();
 478   if(ptD<fMinPtCand) return 0;
 479   if(ptD>fMaxPtCand) return 0;
 480
 481   if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
 482
 483
 484   // selection on daughter tracks
 485   if(selectionLevel==AliRDHFCuts::kAll ||
 486      selectionLevel==AliRDHFCuts::kTracks) {
 487     if(!AreDaughtersSelected(d)) return 0;
 488   }
 489
 490
 491
 492
 493   // selection on candidate
 494   if(selectionLevel==AliRDHFCuts::kAll ||
 495      selectionLevel==AliRDHFCuts::kCandidate) {
 496     //recalculate vertex w/o daughters
 497     AliAODVertex *origownvtx=0x0;
 498     if(fRemoveDaughtersFromPrimary) {
 499       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 500       if(!RecalcOwnPrimaryVtx(d,aod)) {
 501         CleanOwnPrimaryVtx(d,aod,origownvtx);
 502         return 0;
 503       }
 504     }
 505
 506     Int_t okDsKKpi=1;
 507     Int_t okDspiKK=1;
 508     Int_t okMassPhiKKpi=0;
 509     Int_t okMassPhipiKK=0;
 510     Int_t okMassK0starKKpi=0;
 511     Int_t okMassK0starpiKK=0;
 512     Int_t okDsPhiKKpi=0;
 513     Int_t okDsPhipiKK=0;
 514     Int_t okDsK0starKKpi=0;
 515     Int_t okDsK0starpiKK=0;
 516
 517     Double_t pt=d->Pt();
 518     Int_t ptbin=PtBin(pt);
 519     if (ptbin==-1) {
 520       CleanOwnPrimaryVtx(d,aod,origownvtx);
 521       return 0;
 522     }
 523
 524     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
 525     Double_t mDsKKpi=d->InvMassDsKKpi();
 526     Double_t mDspiKK=d->InvMassDspiKK();
 527     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
 528     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
 529     if(!okDsKKpi && !okDspiKK){
 530       CleanOwnPrimaryVtx(d,aod,origownvtx);
 531       return 0;
 532     }
 533
 534
 535
 536     // cuts on resonant decays (via Phi or K0*)
 537     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 538     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 539     if(okDsKKpi){
 540       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
 541       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
 542       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhiKKpi=1;
 543       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starKKpi = 1;
 544       if(!okMassPhiKKpi && !okMassK0starKKpi) okDsKKpi=0;
 545       if(okMassPhiKKpi) okDsPhiKKpi=1;
 546       if(okMassK0starKKpi) okDsK0starKKpi=1;
 547     }
 548     if(okDspiKK){
 549       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
 550       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
 551       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0starpiKK = 1;
 552       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhipiKK=1;
 553       if(!okMassPhipiKK && !okMassK0starpiKK) okDspiKK=0;
 554       if(okMassPhipiKK) okDsPhipiKK=1;
 555       if(okMassK0starpiKK) okDsK0starpiKK=1;
 556     }
 557     if(!okDsKKpi && !okDspiKK){
 558       CleanOwnPrimaryVtx(d,aod,origownvtx);
 559       return 0;
 560     }
 561
 562     // Cuts on track pairs
 563     for(Int_t i=0;i<3;i++){
 564       if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]){
 565         CleanOwnPrimaryVtx(d,aod,origownvtx);
 566         return 0;
 567       }
 568     }
 569     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] ||
 570        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]){
 571       CleanOwnPrimaryVtx(d,aod,origownvtx);
 572       return 0;
 573     }
 574
 575
 576
 577     //single track
 578     if(TMath::Abs(d->Pt2Prong(1)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 579        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]){
 580       CleanOwnPrimaryVtx(d,aod,origownvtx);
 581       return 0;
 582     }
 583
 584     if(okDsKKpi){
 585       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 586          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
 587       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 588          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
 589     }
 590     if(okDspiKK){
 591       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 592          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
 593       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 594          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
 595     }
 596     if(!okDsKKpi && !okDspiKK){
 597       CleanOwnPrimaryVtx(d,aod,origownvtx);
 598       return 0;
 599     }
 600
 601     // Cuts on candidate triplet
 602
 603
 604     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]){
 605       CleanOwnPrimaryVtx(d,aod,origownvtx);
 606       return 0;
 607     }
 608
 609     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 610        d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 611        d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {
 612       CleanOwnPrimaryVtx(d,aod,origownvtx);
 613       return 0;
 614     }
 615
 616     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]){
 617       CleanOwnPrimaryVtx(d,aod,origownvtx);
 618       return 0;
 619     }
 620
 621
 622     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 623     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]){
 624       CleanOwnPrimaryVtx(d,aod,origownvtx);
 625       return 0;
 626     }
 627
 628
 629     //sec vert
 630     Double_t sigmavert=d->GetSigmaVert(aod);
 631     if(sigmavert>fCutsRD[GetGlobalIndex(6,ptbin)]){
 632       CleanOwnPrimaryVtx(d,aod,origownvtx);
 633       return 0;
 634     }
 635
 636     // decay length XY
 637     if(d->DecayLengthXY()<fCutsRD[GetGlobalIndex(16,ptbin)]){
 638       CleanOwnPrimaryVtx(d,aod,origownvtx);
 639       return 0;
 640     }
 641
 642     //norm decay length
 643     if(d->NormalizedDecayLength()<fCutsRD[GetGlobalIndex(17,ptbin)]){
 644       CleanOwnPrimaryVtx(d,aod,origownvtx);
 645       return 0;
 646     }
 647
 648     //norm decay length XY
 649     if(d->NormalizedDecayLengthXY()<fCutsRD[GetGlobalIndex(18,ptbin)]){
 650       CleanOwnPrimaryVtx(d,aod,origownvtx);
 651       return 0;
 652     }
 653
 654     //cos pointing XY
 655     if(d->CosPointingAngleXY()<fCutsRD[GetGlobalIndex(19,ptbin)]){
 656       CleanOwnPrimaryVtx(d,aod,origownvtx);
 657       return 0;
 658     }
 659
 660
 661     if(okDsKKpi){
 662       Double_t cosPiKPhiRFKKpi=d->CosPiKPhiRFrameKKpi();
 663       Double_t kincutPiKPhiKKpi=TMath::Abs(cosPiKPhiRFKKpi*cosPiKPhiRFKKpi*cosPiKPhiRFKKpi);
 664       if(kincutPiKPhiKKpi<fCutsRD[GetGlobalIndex(14,ptbin)]) okDsKKpi=0;
 665     }
 666     if(okDspiKK){
 667       Double_t cosPiKPhiRFpiKK=d->CosPiKPhiRFramepiKK();
 668       Double_t kincutPiKPhipiKK=TMath::Abs(cosPiKPhiRFpiKK*cosPiKPhiRFpiKK*cosPiKPhiRFpiKK);
 669       if(kincutPiKPhipiKK<fCutsRD[GetGlobalIndex(14,ptbin)]) okDspiKK=0;
 670     }
 671     if(!okDsKKpi && !okDspiKK){
 672       CleanOwnPrimaryVtx(d,aod,origownvtx);
 673       return 0;
 674     }
 675
 676
 677
 678     if(okDsKKpi){
 679       Double_t cosPiDsLabFrameKKpi=d->CosPiDsLabFrameKKpi();
 680       if(cosPiDsLabFrameKKpi>fCutsRD[GetGlobalIndex(15,ptbin)]) okDsKKpi=0;
 681     }
 682     if(okDspiKK){
 683       Double_t cosPiDsLabFramepiKK=d->CosPiDsLabFramepiKK();
 684       if(cosPiDsLabFramepiKK>fCutsRD[GetGlobalIndex(15,ptbin)]) okDspiKK=0;
 685     }
 686     if(!okDsKKpi && !okDspiKK){
 687       CleanOwnPrimaryVtx(d,aod,origownvtx);
 688       return 0;
 689     }
 690
 691      // unset recalculated primary vertex when not needed any more
 692     CleanOwnPrimaryVtx(d,aod,origownvtx);
 693
 694
 695
 696     if(!okDsKKpi){
 697       okDsPhiKKpi=0;
 698       okDsK0starKKpi=0;
 699     }
 700     if(!okDspiKK){
 701       okDsPhipiKK=0;
 702       okDsK0starpiKK=0;
 703     }
 704
 705     // PID selection
 706     Int_t returnvaluePID=3;
 707     if(selectionLevel==AliRDHFCuts::kAll ||
 708        selectionLevel==AliRDHFCuts::kCandidate ||
 709        selectionLevel==AliRDHFCuts::kPID) {
 710       returnvaluePID = IsSelectedPID(d);
 711       fIsSelectedPID=returnvaluePID;
 712     }
 713     if(returnvaluePID==0)return 0;
 714
 715     Bool_t okPidDsKKpi=returnvaluePID&1;
 716     Bool_t okPidDspiKK=returnvaluePID&2;
 717     if(!okPidDsKKpi){
 718       okDsPhiKKpi=0;
 719       okDsK0starKKpi=0;
 720     }
 721     if(!okPidDspiKK){
 722       okDsPhipiKK=0;
 723       okDsK0starpiKK=0;
 724     }
 725
 726     if((okPidDsKKpi && okDsKKpi)||(okPidDspiKK && okDspiKK)){
 727       Int_t returnvalue=0;
 728       if(okDsKKpi) returnvalue+=1;
 729       if(okDspiKK) returnvalue+=2;
 730       if(okDsPhiKKpi) returnvalue+=4;
 731       if(okDsPhipiKK) returnvalue+=8;
 732       if(okDsK0starKKpi) returnvalue+=16;
 733       if(okDsK0starpiKK) returnvalue+=32;
 734       return returnvalue;
 735     }else{
 736       return 0;
 737     }
 738   }
 739   return 15;
 740
 741 }
 742
 743 //--------------------------------------------------------------------------
 744
 745 UInt_t AliRDHFCutsDstoKKpi::GetPIDTrackTPCTOFBitMap(AliAODTrack *track) const{
 746
 747   UInt_t bitmap=0;
 748
 749   Double_t sigmaTPCPionHyp=-999.;
 750   Double_t sigmaTPCKaonHyp=-999.;
 751   Double_t sigmaTPCProtonHyp=-999.;
 752   Double_t sigmaTOFPionHyp=-999.;
 753   Double_t sigmaTOFKaonHyp=-999.;
 754   Double_t sigmaTOFProtonHyp=-999.;
 755
 756   Int_t oksigmaTPCPionHyp=fPidHF->GetnSigmaTPC(track,2,sigmaTPCPionHyp);
 757   Int_t oksigmaTPCKaonHyp=fPidHF->GetnSigmaTPC(track,3,sigmaTPCKaonHyp);
 758   Int_t oksigmaTPCProtonHyp=fPidHF->GetnSigmaTPC(track,4,sigmaTPCProtonHyp);
 759   Int_t oksigmaTOFPionHyp=fPidHF->GetnSigmaTOF(track,2,sigmaTOFPionHyp);
 760   Int_t oksigmaTOFKaonHyp=fPidHF->GetnSigmaTOF(track,3,sigmaTOFKaonHyp);
 761   Int_t oksigmaTOFProtonHyp=fPidHF->GetnSigmaTOF(track,4,sigmaTOFProtonHyp);
 762
 763   sigmaTPCPionHyp=TMath::Abs(sigmaTPCPionHyp);
 764   sigmaTPCKaonHyp=TMath::Abs(sigmaTPCKaonHyp);
 765   sigmaTPCProtonHyp=TMath::Abs(sigmaTPCProtonHyp);
 766   sigmaTOFPionHyp=TMath::Abs(sigmaTOFPionHyp);
 767   sigmaTOFKaonHyp=TMath::Abs(sigmaTOFKaonHyp);
 768   sigmaTOFProtonHyp=TMath::Abs(sigmaTOFProtonHyp);
 769
 770  if (oksigmaTPCPionHyp && sigmaTPCPionHyp>0.){
 771     if (sigmaTPCPionHyp<1.) bitmap+=1<<kTPCPionLess1;
 772     else{
 773       if (sigmaTPCPionHyp<2.) bitmap+=1<<kTPCPionMore1Less2;
 774       else {
 775         if (sigmaTPCPionHyp<3.) bitmap+=1<<kTPCPionMore2Less3;
 776         else bitmap+=1<<kTPCPionMore3;
 777       }
 778     }
 779   }
 780
 781   if (oksigmaTPCKaonHyp && sigmaTPCKaonHyp>0.){
 782     if (sigmaTPCKaonHyp<1.) bitmap+=1<<kTPCKaonLess1;
 783     else{
 784       if (sigmaTPCKaonHyp<2.) bitmap+=1<<kTPCKaonMore1Less2;
 785       else {
 786         if (sigmaTPCKaonHyp<3.) bitmap+=1<<kTPCKaonMore2Less3;
 787         else bitmap+=1<<kTPCKaonMore3;
 788       }
 789     }
 790   }
 791
 792   if (oksigmaTPCProtonHyp && sigmaTPCProtonHyp>0.){
 793     if (sigmaTPCProtonHyp<1.) bitmap+=1<<kTPCProtonLess1;
 794     else{
 795       if (sigmaTPCProtonHyp<2.) bitmap+=1<<kTPCProtonMore1Less2;
 796       else {
 797         if (sigmaTPCProtonHyp<3.) bitmap+=1<<kTPCProtonMore2Less3;
 798         else bitmap+=1<<kTPCProtonMore3;
 799       }
 800     }
 801   }
 802
 803   if (oksigmaTOFPionHyp && sigmaTOFPionHyp>0.){
 804     if (sigmaTOFPionHyp<1.) bitmap+=1<<kTOFPionLess1;
 805     else{
 806       if (sigmaTOFPionHyp<2.) bitmap+=1<<kTOFPionMore1Less2;
 807       else {
 808         if (sigmaTOFPionHyp<3.) bitmap+=1<<kTOFPionMore2Less3;
 809         else bitmap+=1<<kTOFPionMore3;
 810       }
 811     }
 812   }
 813
 814   if (oksigmaTOFKaonHyp && sigmaTOFKaonHyp>0.){
 815     if (sigmaTOFKaonHyp<1.) bitmap+=1<<kTOFKaonLess1;
 816     else{
 817       if (sigmaTOFKaonHyp<2.) bitmap+=1<<kTOFKaonMore1Less2;
 818       else {
 819         if (sigmaTOFKaonHyp<3.) bitmap+=1<<kTOFKaonMore2Less3;
 820         else bitmap+=1<<kTOFKaonMore3;
 821       }
 822     }
 823   }
 824
 825   if (oksigmaTOFProtonHyp && sigmaTOFProtonHyp>0.){
 826     if (sigmaTOFProtonHyp<1.) bitmap+=1<<kTOFProtonLess1;
 827     else{
 828       if (sigmaTOFProtonHyp<2.) bitmap+=1<<kTOFProtonMore1Less2;
 829       else {
 830         if (sigmaTOFProtonHyp<3.) bitmap+=1<<kTOFProtonMore2Less3;
 831         else bitmap+=1<<kTOFProtonMore3;
 832       }
 833     }
 834   }
 835
 836
 837
 838   return bitmap;
 839
 840 }
 841
 842
 843 //---------------------------------------------------------------------------
 844
 845 void AliRDHFCutsDstoKKpi::SetStandardCutsPP2010() {
 846   //
 847   //STANDARD CUTS USED FOR 2010 pp analysis
 848   //
 849
 850   SetName("DstoKKpiCutsStandard");
 851   SetTitle("Standard Cuts for D+s analysis");
 852
 853   // PILE UP REJECTION
 854   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 855
 856   // EVENT CUTS
 857   SetMinVtxContr(1);
 858
 859   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
 860   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 861   //default
 862   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 863   esdTrackCuts->SetRequireITSRefit(kTRUE);
 864   //esdTrackCuts->SetMinNClustersITS(4); // default is 5
 865   esdTrackCuts->SetMinNClustersTPC(70);
 866   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 867                                          AliESDtrackCuts::kAny);
 868   // default is kBoth, otherwise kAny
 869   esdTrackCuts->SetMinDCAToVertexXY(0.);
 870   esdTrackCuts->SetPtRange(0.3,1.e10);
 871
 872   AddTrackCuts(esdTrackCuts);
 873
 874
 875
 876   const Int_t nptbins=4;
 877   Float_t ptbins[nptbins+1];
 878   ptbins[0]=2.;
 879   ptbins[1]=4.;
 880   ptbins[2]=6.;
 881   ptbins[3]=8.;
 882   ptbins[4]=12.;
 883
 884   const Int_t nvars=20;
 885
 886   Float_t** anacutsval;
 887   anacutsval=new Float_t*[nvars];
 888
 889   for(Int_t ic=0;ic<nvars;ic++){anacutsval[ic]=new Float_t[nptbins];}
 890   for(Int_t ipt=0;ipt<nptbins;ipt++){
 891
 892     anacutsval[0][ipt]=0.35;
 893     anacutsval[1][ipt]=0.3;
 894     anacutsval[2][ipt]=0.3;
 895     anacutsval[3][ipt]=0.;
 896     anacutsval[4][ipt]=0.;
 897     anacutsval[5][ipt]=0.005;
 898     anacutsval[8][ipt]=0.;
 899     anacutsval[10][ipt]=0.;
 900     anacutsval[11][ipt]=1000.0;
 901     anacutsval[13][ipt]=0.1;
 902     anacutsval[16][ipt]=0.;
 903     anacutsval[17][ipt]=0.;
 904     anacutsval[18][ipt]=0.;
 905     anacutsval[19][ipt]=-1.;
 906
 907
 908   }
 909
 910   //sigmavertex
 911
 912   anacutsval[6][0]=0.020;
 913   anacutsval[6][1]=0.030;
 914   anacutsval[6][2]=0.030;
 915   anacutsval[6][3]=0.060;
 916
 917   //Decay length
 918
 919   anacutsval[7][0]=0.035;
 920   anacutsval[7][1]=0.035;
 921   anacutsval[7][2]=0.040;
 922   anacutsval[7][3]=0.040;
 923
 924   //cosThetaPoint
 925
 926   anacutsval[9][0]=0.94;
 927   anacutsval[9][1]=0.94;
 928   anacutsval[9][2]=0.94;
 929   anacutsval[9][3]=0.94;
 930
 931   //phi DeltaMass
 932
 933   anacutsval[12][0]=0.0080;
 934   anacutsval[12][1]=0.0050;
 935   anacutsval[12][2]=0.0045;
 936   anacutsval[12][3]=0.0090;
 937
 938   //Kin1
 939
 940   anacutsval[14][0]=0.10;
 941   anacutsval[14][1]=0.05;
 942   anacutsval[14][2]=0.0;
 943   anacutsval[14][3]=0.05;
 944
 945   //Kin2
 946
 947   anacutsval[15][0]=0.95;
 948   anacutsval[15][1]=0.95;
 949   anacutsval[15][2]=1.;
 950   anacutsval[15][3]=0.95;
 951
 952   fPidHF->SetOldPid(kTRUE);
 953   SetUsePID(kTRUE);
 954   SetPidOption(1);
 955   SetMaxPtStrongPid(9999.);
 956   SetGlobalIndex(nvars,nptbins);
 957   SetPtBins(nptbins+1,ptbins);
 958   SetCuts(nvars,nptbins,anacutsval);
 959   SetRemoveDaughtersFromPrim(kTRUE);
 960
 961   PrintAll();
 962
 963   for(Int_t iic=0;iic<nvars;iic++){delete [] anacutsval[iic];}
 964   delete [] anacutsval;
 965   anacutsval=NULL;
 966
 967   return;
 968 }
 969