PWG3/vertexingHF/AliRDHFCutsD0toKpi.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 D0->Kpi
  21 //
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TDatabasePDG.h>
  26 #include <Riostream.h>
  27
  28 #include "AliRDHFCutsD0toKpi.h"
  29 #include "AliAODRecoDecayHF2Prong.h"
  30 #include "AliAODTrack.h"
  31 #include "AliESDtrack.h"
  32 #include "AliAODPid.h"
  33 #include "AliTPCPIDResponse.h"
  34 #include "AliAODVertex.h"
  35 #include "AliKFVertex.h"
  36 #include "AliKFParticle.h"
  37
  38 ClassImp(AliRDHFCutsD0toKpi)
  39
  40 //--------------------------------------------------------------------------
  41 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) :
  42 AliRDHFCuts(name),
  43 fUseSpecialCuts(kFALSE),
  44 fLowPt(kTRUE),
  45 fDefaultPID(kTRUE),
  46 fUseKF(kFALSE)
  47 {
  48   //
  49   // Default Constructor
  50   //
  51   Int_t nvars=9;
  52   SetNVars(nvars);
  53   TString varNames[9]={"inv. mass [GeV]",
  54                        "dca [cm]",
  55                        "cosThetaStar",
  56                        "pTK [GeV/c]",
  57                        "pTPi [GeV/c]",
  58                        "d0K [cm]",
  59                        "d0Pi [cm]",
  60                        "d0d0 [cm^2]",
  61                        "cosThetaPoint"};
  62   Bool_t isUpperCut[9]={kTRUE,
  63                         kTRUE,
  64                         kTRUE,
  65                         kFALSE,
  66                         kFALSE,
  67                         kTRUE,
  68                         kTRUE,
  69                         kTRUE,
  70                         kFALSE};
  71   SetVarNames(nvars,varNames,isUpperCut);
  72   Bool_t forOpt[9]={kFALSE,
  73                     kTRUE,
  74                     kTRUE,
  75                     kFALSE,
  76                     kFALSE,
  77                     kFALSE,
  78                     kFALSE,
  79                     kTRUE,
  80                     kTRUE};
  81   SetVarsForOpt(4,forOpt);
  82   Float_t limits[2]={0,999999999.};
  83   SetPtBins(2,limits);
  84
  85 }
  86 //--------------------------------------------------------------------------
  87 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
  88   AliRDHFCuts(source),
  89   fUseSpecialCuts(source.fUseSpecialCuts),
  90   fLowPt(source.fLowPt),
  91   fDefaultPID(source.fDefaultPID),
  92   fUseKF(source.fUseKF)
  93 {
  94   //
  95   // Copy constructor
  96   //
  97
  98 }
  99 //--------------------------------------------------------------------------
 100 AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &source)
 101 {
 102   //
 103   // assignment operator
 104   //
 105   if(&source == this) return *this;
 106
 107   AliRDHFCuts::operator=(source);
 108   fUseSpecialCuts=source.fUseSpecialCuts;
 109   fLowPt=source.fLowPt;
 110   fDefaultPID=source.fDefaultPID;
 111
 112   return *this;
 113 }
 114
 115
 116 //---------------------------------------------------------------------------
 117 void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
 118   //
 119   // Fills in vars the values of the variables
 120   //
 121
 122   if(nvars!=fnVarsForOpt) {
 123     printf("AliRDHFCutsD0toKpi::GetCutsVarsForOpt: wrong number of variables\n");
 124     return;
 125   }
 126
 127   AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
 128
 129   Int_t iter=-1;
 130   if(fVarsForOpt[0]){
 131     iter++;
 132     if(TMath::Abs(pdgdaughters[0])==211) {
 133       vars[iter]=dd->InvMassD0();
 134     } else {
 135       vars[iter]=dd->InvMassD0bar();
 136     }
 137   }
 138   if(fVarsForOpt[1]){
 139     iter++;
 140     vars[iter]=dd->GetDCA();
 141   }
 142   if(fVarsForOpt[2]){
 143     iter++;
 144     if(TMath::Abs(pdgdaughters[0])==211) {
 145       vars[iter] = dd->CosThetaStarD0();
 146     } else {
 147       vars[iter] = dd->CosThetaStarD0bar();
 148     }
 149   }
 150   if(fVarsForOpt[3]){
 151     iter++;
 152    if(TMath::Abs(pdgdaughters[0])==321) {
 153      vars[iter]=dd->PtProng(0);
 154    }
 155    else{
 156      vars[iter]=dd->PtProng(1);
 157    }
 158   }
 159   if(fVarsForOpt[4]){
 160     iter++;
 161    if(TMath::Abs(pdgdaughters[0])==211) {
 162      vars[iter]=dd->PtProng(0);
 163    }
 164    else{
 165      vars[iter]=dd->PtProng(1);
 166    }
 167   }
 168   if(fVarsForOpt[5]){
 169     iter++;
 170     if(TMath::Abs(pdgdaughters[0])==321) {
 171      vars[iter]=dd->Getd0Prong(0);
 172    }
 173    else{
 174      vars[iter]=dd->Getd0Prong(1);
 175    }
 176   }
 177   if(fVarsForOpt[6]){
 178     iter++;
 179      if(TMath::Abs(pdgdaughters[0])==211) {
 180      vars[iter]=dd->Getd0Prong(0);
 181    }
 182    else{
 183      vars[iter]=dd->Getd0Prong(1);
 184    }
 185   }
 186   if(fVarsForOpt[7]){
 187     iter++;
 188     vars[iter]= dd->Prodd0d0();
 189   }
 190   if(fVarsForOpt[8]){
 191     iter++;
 192     vars[iter]=dd->CosPointingAngle();
 193   }
 194
 195   return;
 196 }
 197 //---------------------------------------------------------------------------
 198 Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent* aod) {
 199   //
 200   // Apply selection
 201   //
 202
 203
 204   fIsSelectedCuts=0;
 205   fIsSelectedPID=0;
 206
 207   if(!fCutsRD){
 208     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 209     return 0;
 210   }
 211   //PrintAll();
 212   AliAODRecoDecayHF2Prong* d=(AliAODRecoDecayHF2Prong*)obj;
 213
 214   if(!d){
 215     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 216     return 0;
 217   }
 218
 219   Double_t ptD=d->Pt();
 220   if(ptD<fMinPtCand) return 0;
 221   if(ptD>fMaxPtCand) return 0;
 222
 223   // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
 224   Int_t returnvaluePID=3;
 225   Int_t returnvalueCuts=3;
 226
 227   // selection on candidate
 228   if(selectionLevel==AliRDHFCuts::kAll ||
 229      selectionLevel==AliRDHFCuts::kCandidate) {
 230
 231     if(!fUseKF) {
 232
 233       //recalculate vertex w/o daughters
 234       AliAODVertex *origownvtx=0x0;
 235       if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
 236         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 237         if(!RecalcOwnPrimaryVtx(d,aod)) {
 238           CleanOwnPrimaryVtx(d,aod,origownvtx);
 239           return 0;
 240         }
 241       }
 242
 243       if(fUseMCVertex) {
 244         if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 245         if(!SetMCPrimaryVtx(d,aod)) {
 246           CleanOwnPrimaryVtx(d,aod,origownvtx);
 247           return 0;
 248         }
 249       }
 250
 251       Double_t pt=d->Pt();
 252
 253       Int_t okD0=0,okD0bar=0;
 254
 255       Int_t ptbin=PtBin(pt);
 256       if (ptbin==-1) {
 257         CleanOwnPrimaryVtx(d,aod,origownvtx);
 258         return 0;
 259       }
 260
 261       Double_t mD0,mD0bar,ctsD0,ctsD0bar;
 262       okD0=1; okD0bar=1;
 263
 264       Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 265
 266       d->InvMassD0(mD0,mD0bar);
 267       if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
 268       if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)])  okD0bar = 0;
 269       if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 270
 271       if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 272
 273
 274       if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
 275       if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
 276       if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 277
 278
 279       if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
 280          TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
 281       if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
 282          TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
 283       if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 284
 285       if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 286
 287
 288       d->CosThetaStarD0(ctsD0,ctsD0bar);
 289       if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
 290       if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
 291       if(!okD0 && !okD0bar)   {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 292
 293       if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 294
 295       if (returnvalueCuts!=0) {
 296         if (okD0) returnvalueCuts=1; //cuts passed as D0
 297         if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
 298         if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
 299       }
 300
 301       // call special cuts
 302       Int_t special=1;
 303       if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
 304       //if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
 305
 306       // unset recalculated primary vertex when not needed any more
 307       CleanOwnPrimaryVtx(d,aod,origownvtx);
 308
 309     } else {
 310       // go to selection with Kalman vertexing, if requested
 311       returnvalueCuts = IsSelectedKF(d,aod);
 312     }
 313
 314     fIsSelectedCuts=returnvalueCuts;
 315     if(!returnvalueCuts) return 0;
 316   }
 317
 318
 319   // selection on PID
 320   if(selectionLevel==AliRDHFCuts::kAll ||
 321      selectionLevel==AliRDHFCuts::kCandidate ||
 322      selectionLevel==AliRDHFCuts::kPID) {
 323     returnvaluePID = IsSelectedPID(d);
 324     fIsSelectedPID=returnvaluePID;
 325     if(!returnvaluePID) return 0;
 326   }
 327
 328   Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
 329
 330   if(!returnvalueComb) return 0;
 331
 332   // selection on daughter tracks
 333   if(selectionLevel==AliRDHFCuts::kAll ||
 334      selectionLevel==AliRDHFCuts::kTracks) {
 335     if(!AreDaughtersSelected(d)) return 0;
 336   }
 337
 338   //  cout<<"Pid = "<<returnvaluePID<<endl;
 339   return returnvalueComb;
 340 }
 341
 342 //------------------------------------------------------------------------------------------
 343 Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
 344                                        AliAODEvent* aod) const {
 345   //
 346   // Apply selection using KF-vertexing
 347   //
 348
 349   AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
 350   AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1);
 351
 352   if(!track0 || !track1) {
 353     cout<<"one or two D0 daughters missing!"<<endl;
 354     return 0;
 355   }
 356
 357   // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
 358   Int_t returnvalueCuts=3;
 359
 360   // candidate selection with AliKF
 361   AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
 362
 363   Int_t okD0=0,okD0bar=0;
 364   okD0=1; okD0bar=1;
 365
 366   // convert tracks into AliKFParticles
 367
 368   AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
 369   AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
 370   AliKFParticle posPiKF(*track0,211); // pos pion kandidate
 371   AliKFParticle posKKF(*track0,321); // pos kaon kandidate
 372
 373   // build D0 candidates
 374
 375   AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
 376   AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
 377
 378   // create kf primary vertices
 379
 380   AliAODVertex *vtx1 = aod->GetPrimaryVertex();
 381   AliKFVertex primVtx(*vtx1);
 382   AliKFVertex aprimVtx(*vtx1);
 383
 384   if(primVtx.GetNContributors()<=0) okD0 = 0;
 385   if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
 386   if(!okD0 && !okD0bar) returnvalueCuts=0;
 387
 388   // calculate mass
 389
 390   Double_t d0mass = d0c.GetMass();
 391   Double_t ad0mass = ad0c.GetMass();
 392
 393   // calculate P of D0 and D0bar
 394   Double_t d0P = d0c.GetP();
 395   Double_t d0Px = d0c.GetPx();
 396   Double_t d0Py = d0c.GetPy();
 397   Double_t d0Pz = d0c.GetPz();
 398   Double_t ad0P = ad0c.GetP();
 399   Double_t ad0Px = ad0c.GetPx();
 400   Double_t ad0Py = ad0c.GetPy();
 401   Double_t ad0Pz = ad0c.GetPz();
 402
 403   //calculate Pt of D0 and D0bar
 404
 405   Double_t pt=d0c.GetPt();
 406   Double_t apt=ad0c.GetPt();
 407
 408   // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
 409
 410   if(track0->GetUsedForPrimVtxFit()) {
 411     primVtx -= posPiKF;
 412     aprimVtx -= posKKF;
 413   }
 414
 415   if(track1->GetUsedForPrimVtxFit()) {
 416     primVtx -= negKKF;
 417     aprimVtx -= negPiKF;
 418   }
 419
 420   primVtx += d0c;
 421   aprimVtx += ad0c;
 422
 423   if(primVtx.GetNContributors()<=0) okD0 = 0;
 424   if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
 425   if(!okD0 && !okD0bar) returnvalueCuts=0;
 426
 427   //calculate cut variables
 428
 429   // calculate impact params of daughters w.r.t recalculated vertices
 430
 431   Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
 432   Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
 433   Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
 434   Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
 435
 436   // calculate Product of Impact Params
 437
 438   Double_t prodParam = impactPi*impactKa;
 439   Double_t aprodParam = aimpactPi*aimpactKa;
 440
 441   // calculate cosine of pointing angles
 442
 443   TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
 444   TVector3 fline(d0c.GetX()-primVtx.GetX(),
 445                  d0c.GetY()-primVtx.GetY(),
 446                  d0c.GetZ()-primVtx.GetZ());
 447   Double_t pta = mom.Angle(fline);
 448   Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
 449
 450   TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
 451   TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
 452                   ad0c.GetY()-aprimVtx.GetY(),
 453                   ad0c.GetZ()-aprimVtx.GetZ());
 454   Double_t apta = amom.Angle(afline);
 455   Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
 456
 457   // calculate P of Pions at Decay Position of D0 and D0bar candidates
 458   negKKF.TransportToParticle(d0c);
 459   posPiKF.TransportToParticle(d0c);
 460   posKKF.TransportToParticle(ad0c);
 461   negPiKF.TransportToParticle(ad0c);
 462
 463   Double_t pxPi =  posPiKF.GetPx();
 464   Double_t pyPi =  posPiKF.GetPy();
 465   Double_t pzPi =  posPiKF.GetPz();
 466   Double_t ptPi =  posPiKF.GetPt();
 467
 468   Double_t apxPi =  negPiKF.GetPx();
 469   Double_t apyPi =  negPiKF.GetPy();
 470   Double_t apzPi =  negPiKF.GetPz();
 471   Double_t aptPi =  negPiKF.GetPt();
 472
 473   // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
 474
 475   Double_t ptK =  negKKF.GetPt();
 476   Double_t aptK =  posKKF.GetPt();
 477
 478   //calculate cos(thetastar)
 479   Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 480   Double_t massp[2];
 481   massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
 482   massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
 483   Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
 484                                -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
 485
 486   // cos(thetastar) for D0 and Pion
 487
 488   Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
 489   Double_t beta = d0P/d0E;
 490   Double_t gamma = d0E/massvtx;
 491   TVector3 momPi(pxPi,pyPi,pzPi);
 492   TVector3 momTot(d0Px,d0Py,d0Pz);
 493   Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
 494   Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
 495
 496   // cos(thetastar) for D0bar and Pion
 497
 498   Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
 499   Double_t abeta = ad0P/ad0E;
 500   Double_t agamma = ad0E/massvtx;
 501   TVector3 amomPi(apxPi,apyPi,apzPi);
 502   TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
 503   Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
 504   Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;
 505
 506   // calculate reduced Chi2 for the full D0 fit
 507   d0c.SetProductionVertex(primVtx);
 508   ad0c.SetProductionVertex(aprimVtx);
 509   negKKF.SetProductionVertex(d0c);
 510   posPiKF.SetProductionVertex(d0c);
 511   posKKF.SetProductionVertex(ad0c);
 512   negPiKF.SetProductionVertex(ad0c);
 513   d0c.TransportToProductionVertex();
 514   ad0c.TransportToProductionVertex();
 515
 516   // calculate the decay length
 517   Double_t decayLengthD0 = d0c.GetDecayLength();
 518   Double_t adecayLengthD0 = ad0c.GetDecayLength();
 519
 520   Double_t chi2D0 = 50.;
 521   if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
 522     chi2D0 = d0c.GetChi2()/d0c.GetNDF();
 523   }
 524
 525   Double_t achi2D0 = 50.;
 526   if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
 527     achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
 528   }
 529
 530   // Get the Pt-bins
 531   Int_t ptbin=PtBin(pt);
 532   Int_t aptbin=PtBin(apt);
 533
 534   if(ptbin < 0) okD0 = 0;
 535   if(aptbin < 0) okD0bar = 0;
 536   if(!okD0 && !okD0bar) returnvalueCuts=0;
 537
 538   if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
 539   if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
 540   if(!okD0 && !okD0bar) returnvalueCuts=0;
 541
 542
 543   if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
 544      TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
 545
 546   if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
 547      TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
 548
 549   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 550
 551   // for the moment via the standard method due to bug in AliKF
 552   if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
 553   if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
 554   if(!okD0 && !okD0bar) returnvalueCuts=0;
 555
 556
 557   if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
 558   if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)])  okD0bar = 0;
 559   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 560
 561
 562   if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
 563   if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
 564   if(!okD0 && !okD0bar)   returnvalueCuts=0;
 565
 566   if(prodParam  > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
 567   if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
 568   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 569
 570   if(cosP  < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0;
 571   if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
 572   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 573
 574   if(chi2D0  > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0;
 575   if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
 576   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 577
 578   if(decayLengthD0  < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0;
 579   if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
 580   if(!okD0 && !okD0bar)  returnvalueCuts=0;
 581
 582   if(returnvalueCuts!=0) {
 583     if(okD0) returnvalueCuts=1; //cuts passed as D0
 584     if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
 585     if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
 586   }
 587
 588   return returnvalueCuts;
 589 }
 590
 591 //---------------------------------------------------------------------------
 592
 593 Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 594 {
 595   //
 596   // Checking if D0 is in fiducial acceptance region
 597   //
 598
 599   if(pt > 5.) {
 600     // applying cut for pt > 5 GeV
 601     AliDebug(2,Form("pt of D0 = %f (> 5), cutting at |y| < 0.8\n",pt));
 602     if (TMath::Abs(y) > 0.8){
 603       return kFALSE;
 604     }
 605   } else {
 606     // appliying smooth cut for pt < 5 GeV
 607     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 608     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 609     AliDebug(2,Form("pt of D0 = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 610     if (y < minFiducialY || y > maxFiducialY){
 611       return kFALSE;
 612     }
 613   }
 614
 615   return kTRUE;
 616 }
 617 //---------------------------------------------------------------------------
 618 Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
 619 {
 620   // ############################################################
 621   //
 622   // Apply PID selection
 623   //
 624   //
 625   // ############################################################
 626
 627   if(!fUsePID) return 3;
 628   if(fDefaultPID) return IsSelectedPIDdefault(d);
 629   fWhyRejection=0;
 630   Int_t isD0D0barPID[2]={1,2};
 631   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
 632   //                                                                                                 same for prong 2
 633   //                                               values convention -1 = discarded
 634   //                                                                  0 = not identified (but compatible) || No PID (->hasPID flag)
 635   //                                                                  1 = identified
 636   // PID search:   pion (TPC) or not K (TOF), Kaon hypothesis for both
 637   // Initial hypothesis: unknwon (but compatible)
 638   combinedPID[0][0]=0;  // prima figlia, Kaon
 639   combinedPID[0][1]=0;  // prima figlia, pione
 640   combinedPID[1][0]=0;  // seconda figlia, Kaon
 641   combinedPID[1][1]=0;  // seconda figlia, pion
 642
 643   Bool_t checkPIDInfo[2]={kTRUE,kTRUE};
 644   Double_t sigma_tmp[3]={fPidHF->GetSigma(0),fPidHF->GetSigma(1),fPidHF->GetSigma(2)};
 645   for(Int_t daught=0;daught<2;daught++){
 646     //Loop con prongs
 647     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
 648     if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
 649
 650     if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
 651       checkPIDInfo[daught]=kFALSE;
 652       continue;
 653     }
 654
 655     // identify kaon
 656     combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 657
 658     // identify pion
 659
 660     if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
 661      combinedPID[daught][1]=0;
 662     }else{
 663       fPidHF->SetTOF(kFALSE);
 664       combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
 665       fPidHF->SetTOF(kTRUE);
 666       fPidHF->SetCompat(kTRUE);
 667      }
 668
 669
 670    if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
 671     isD0D0barPID[0]=0;
 672     isD0D0barPID[1]=0;
 673    }
 674    else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
 675     if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
 676     else isD0D0barPID[0]=0;// if K+ D0 excluded
 677    }
 678    else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
 679     isD0D0barPID[0]=0;
 680     isD0D0barPID[1]=0;
 681    }
 682    else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
 683    if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
 684    else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
 685         }
 686    else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
 687     if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
 688     else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
 689    }
 690
 691     if(fLowPt && d->Pt()<2.){
 692      Double_t sigmaTPC[3]={3.,2.,0.};
 693      fPidHF->SetSigmaForTPC(sigmaTPC);
 694     // identify kaon
 695     combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 696
 697     Double_t ptProng=aodtrack->P();
 698
 699     if(ptProng<0.6){
 700      fPidHF->SetCompat(kFALSE);
 701      combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 702      fPidHF->SetCompat(kTRUE);
 703     }
 704
 705     if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
 706      combinedPID[daught][1]=0;
 707     }else{
 708       fPidHF->SetTOF(kFALSE);
 709       Double_t sigmaTPCpi[3]={3.,3.,0.};
 710       fPidHF->SetSigmaForTPC(sigmaTPCpi);
 711       combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
 712       fPidHF->SetTOF(kTRUE);
 713        if(ptProng<0.8){
 714         Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
 715         if(isTPCpion){
 716          combinedPID[daught][1]=1;
 717         }else{
 718          combinedPID[daught][1]=-1;
 719         }
 720       }
 721     }
 722
 723    }
 724    fPidHF->SetSigmaForTPC(sigma_tmp);
 725   }// END OF LOOP ON DAUGHTERS
 726
 727    if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
 728     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 729     return 0;
 730    }
 731
 732
 733   // FURTHER PID REQUEST (both daughter info is needed)
 734   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
 735     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
 736     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 737     return 0;
 738   }
 739
 740   if(d->Pt()<2.){
 741     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 742     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
 743       fWhyRejection=32;// reject cases where the Kaon is not identified
 744       return 0;
 745     }
 746   }
 747     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 748
 749   //  cout<<"Why? "<<fWhyRejection<<endl;
 750   return isD0D0barPID[0]+isD0D0barPID[1];
 751 }
 752 //---------------------------------------------------------------------------
 753 Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
 754 {
 755   // ############################################################
 756   //
 757   // Apply PID selection
 758   //
 759   //
 760   // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
 761   //
 762   // d must be a AliAODRecoDecayHF2Prong object
 763   // returns 0 if both D0 and D0bar are rejectecd
 764   //         1 if D0 is accepted while D0bar is rejected
 765   //         2 if D0bar is accepted while D0 is rejected
 766   //         3 if both are accepted
 767   // fWhyRejection variable (not returned for the moment, print it if needed)
 768   //               keeps some information on why a candidate has been
 769   //               rejected according to the following (unfriendly?) scheme
 770   //             if more rejection cases are considered interesting, just add numbers
 771   //
 772   //      TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
 773   //              from 20 to 30: "detector" selection (PID acceptance)
 774   //                                                  26: TPC refit
 775   //                                                  27: ITS refit
 776   //                                                  28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
 777   //
 778   //              from 30 to 40: PID selection
 779   //                                                  31: no Kaon compatible tracks found between daughters
 780   //                                                  32: no Kaon identified tracks found (strong sel. at low momenta)
 781   //                                                  33: both mass hypotheses are rejected
 782   //
 783   // ############################################################
 784
 785   if(!fUsePID) return 3;
 786   fWhyRejection=0;
 787   Int_t isD0D0barPID[2]={1,2};
 788   Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
 789   Double_t tofSig,times[5];// used fot TOF pid
 790   Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
 791   Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
 792   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
 793   //                                                                                                 same for prong 2
 794   //                                               values convention -1 = discarded
 795   //                                                                  0 = not identified (but compatible) || No PID (->hasPID flag)
 796   //                                                                  1 = identified
 797   // PID search:   pion (TPC) or not K (TOF), Kaon hypothesis for both
 798   // Initial hypothesis: unknwon (but compatible)
 799   isKaonPionTOF[0][0]=0;
 800   isKaonPionTOF[0][1]=0;
 801   isKaonPionTOF[1][0]=0;
 802   isKaonPionTOF[1][1]=0;
 803
 804   isKaonPionTPC[0][0]=0;
 805   isKaonPionTPC[0][1]=0;
 806   isKaonPionTPC[1][0]=0;
 807   isKaonPionTPC[1][1]=0;
 808
 809   combinedPID[0][0]=0;
 810   combinedPID[0][1]=0;
 811   combinedPID[1][0]=0;
 812   combinedPID[1][1]=0;
 813
 814
 815
 816   for(Int_t daught=0;daught<2;daught++){
 817     //Loop con prongs
 818
 819     // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
 820
 821     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
 822
 823     if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
 824       fWhyRejection=26;
 825       return 0;
 826     }
 827     if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
 828       fWhyRejection=27;
 829       return 0;
 830     }
 831
 832     AliAODPid *pid=aodtrack->GetDetPid();
 833     if(!pid) {
 834       //delete esdtrack;
 835       hasPID[daught]--;
 836       continue;
 837     }
 838
 839     // ########### Step 1- Check of TPC and TOF response ####################
 840
 841     Double_t ptrack=aodtrack->P();
 842     //#################### TPC PID #######################
 843      if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
 844        // NO TPC PID INFO FOR THIS TRACK
 845        hasPID[daught]--;
 846      }
 847      else {
 848        static AliTPCPIDResponse theTPCpid;
 849        AliAODPid *pidObj = aodtrack->GetDetPid();
 850        Double_t ptProng=pidObj->GetTPCmomentum();
 851        nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
 852        nsigmaTPCK =  theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
 853        //if(ptrack<0.6){
 854        if(ptProng<0.6){
 855          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
 856          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 857          if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
 858          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 859        }
 860        //else if(ptrack<.8){
 861        else if(ptProng<.8){
 862          if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
 863          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 864          if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
 865          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 866        }
 867        else {
 868          //     if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
 869          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 870          //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
 871          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 872        }
 873      }
 874
 875
 876     // ##### TOF PID: do not ask nothing for pion/protons ############
 877      if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
 878        // NO TOF PID INFO FOR THIS TRACK
 879        hasPID[daught]--;
 880      }
 881      else{
 882        tofSig=pid->GetTOFsignal();
 883        pid->GetIntegratedTimes(times);
 884        if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
 885        if(TMath::Abs(tofSig-times[3])>3.*160.){
 886          isKaonPionTOF[daught][0]=-1;
 887        }
 888        else {
 889          if(ptrack<1.5){
 890            isKaonPionTOF[daught][0]=1;
 891          }
 892        }
 893      }
 894
 895      //######### Step 2: COMBINE TOF and TPC PID ###############
 896      // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
 897      combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
 898      combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
 899
 900
 901      //######### Step 3:   USE PID INFO
 902
 903      if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
 904        isD0D0barPID[0]=0;
 905        isD0D0barPID[1]=0;
 906      }
 907      else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-), if k for both TPC and TOF -> is K
 908        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
 909        else isD0D0barPID[0]=0;// if K+ D0 excluded
 910      }
 911      else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){// if in conflict (both pi- and K-) and k- only for TPC or TOF -> reject
 912        isD0D0barPID[0]=0;
 913        isD0D0barPID[1]=0;
 914      }
 915      else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
 916        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
 917        else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
 918      }
 919      else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
 920        if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
 921       else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
 922      }
 923
 924      // ##########  ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification      ###############################
 925      // ########## more tolerant criteria for single particle ID-> more selective criteria for D0   ##############################
 926      // ###############                     NOT OPTIMIZED YET                                  ###################################
 927      if(d->Pt()<2.){
 928        isKaonPionTPC[daught][0]=0;
 929        isKaonPionTPC[daught][1]=0;
 930        AliAODPid *pidObj = aodtrack->GetDetPid();
 931        Double_t ptProng=pidObj->GetTPCmomentum();
 932        //if(ptrack<0.6){
 933        if(ptProng<0.6){
 934          if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
 935          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 936          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
 937          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 938      }
 939        //else if(ptrack<.8){
 940        else if(ptProng<.8){
 941          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
 942          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 943          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
 944          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 945        }
 946        else {
 947          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 948          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 949        }
 950      }
 951
 952   }// END OF LOOP ON DAUGHTERS
 953
 954   // FURTHER PID REQUEST (both daughter info is needed)
 955   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
 956     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
 957     return 0;
 958   }
 959   else if(hasPID[0]==0&&hasPID[1]==0){
 960     fWhyRejection=28;// reject cases in which no PID info is available
 961     return 0;
 962   }
 963   if(d->Pt()<2.){
 964     // request of K identification at low D0 pt
 965     combinedPID[0][0]=0;
 966     combinedPID[0][1]=0;
 967     combinedPID[1][0]=0;
 968     combinedPID[1][1]=0;
 969
 970     combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
 971     combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
 972     combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
 973     combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
 974
 975     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
 976       fWhyRejection=32;// reject cases where the Kaon is not identified
 977       return 0;
 978     }
 979   }
 980
 981   //  cout<<"Why? "<<fWhyRejection<<endl;
 982   return isD0D0barPID[0]+isD0D0barPID[1];
 983 }
 984
 985
 986
 987 //---------------------------------------------------------------------------
 988 Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
 989                                                  Int_t selectionvalCand,
 990                                                  Int_t selectionvalPID) const
 991 {
 992   //
 993   // This method combines the tracks, PID and cuts selection results
 994   //
 995   if(selectionvalTrack==0) return 0;
 996
 997   Int_t returnvalue;
 998
 999   switch(selectionvalPID) {
1000   case 0:
1001     returnvalue=0;
1002     break;
1003   case 1:
1004     returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1005     break;
1006   case 2:
1007     returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1008     break;
1009   case 3:
1010     returnvalue=selectionvalCand;
1011     break;
1012   default:
1013     returnvalue=0;
1014     break;
1015   }
1016
1017   return returnvalue;
1018 }
1019 //----------------------------------------------------------------------------
1020 Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1021 {
1022   //
1023   // Note: this method is temporary
1024   // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1025   //
1026
1027   //apply cuts
1028
1029   Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1030   // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1031   // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1032
1033   Int_t returnvalue=3; //cut passed
1034   for(Int_t i=0;i<2/*prongs*/;i++){
1035     if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1036   }
1037   if(d->DecayLength2()<decLengthCut*decLengthCut)  return 0; //decLengthCut not passed
1038   if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut)  return 0; //decLengthCut not passed
1039
1040
1041   return returnvalue;
1042 }
1043
1044 //----------------------------------------------
1045 void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1046 {
1047   //switch on candidate selection via AliKFparticle
1048   if(!useKF) return;
1049   if(useKF){
1050     fUseKF=useKF;
1051     Int_t nvarsKF=11;
1052     SetNVars(nvarsKF);
1053     TString varNamesKF[11]={"inv. mass [GeV]",
1054                             "dca [cm]",
1055                             "cosThetaStar",
1056                             "pTK [GeV/c]",
1057                             "pTPi [GeV/c]",
1058                             "d0K [cm]",
1059                             "d0Pi [cm]",
1060                             "d0d0 [cm^2]",
1061                             "cosThetaPoint"
1062                             "DecayLength[cm]",
1063                             "RedChi2"};
1064     Bool_t isUpperCutKF[11]={kTRUE,
1065                              kTRUE,
1066                              kTRUE,
1067                              kFALSE,
1068                              kFALSE,
1069                              kTRUE,
1070                              kTRUE,
1071                              kTRUE,
1072                              kFALSE,
1073                              kFALSE,
1074                              kTRUE};
1075     SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1076     SetGlobalIndex();
1077     Bool_t forOpt[11]={kFALSE,
1078                     kTRUE,
1079                     kTRUE,
1080                     kFALSE,
1081                     kFALSE,
1082                     kFALSE,
1083                     kFALSE,
1084                     kTRUE,
1085                     kTRUE,
1086                     kFALSE,
1087                     kFALSE};
1088     SetVarsForOpt(4,forOpt);
1089   }
1090   return;
1091 }
1092
1093
1094 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1095   //
1096   //STANDARD CUTS USED FOR 2010 pp analysis
1097   //dca cut will be enlarged soon to 400 micron
1098   //
1099
1100   SetName("D0toKpiCutsStandard");
1101   SetTitle("Standard Cuts for D0 analysis");
1102
1103   // PILE UP REJECTION
1104   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1105
1106   // EVENT CUTS
1107   SetMinVtxContr(1);
1108
1109
1110   // TRACKS ON SINGLE TRACKS
1111   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1112   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1113   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1114   esdTrackCuts->SetRequireITSRefit(kTRUE);
1115   //  esdTrackCuts->SetMinNClustersITS(4);
1116   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1117   esdTrackCuts->SetMinDCAToVertexXY(0.);
1118   esdTrackCuts->SetEtaRange(-0.8,0.8);
1119   esdTrackCuts->SetPtRange(0.3,1.e10);
1120
1121   AddTrackCuts(esdTrackCuts);
1122
1123   const Int_t nptbins =13;
1124   const Double_t ptmax = 9999.;
1125   const Int_t nvars=9;
1126   Float_t ptbins[nptbins+1];
1127   ptbins[0]=0.;
1128   ptbins[1]=0.5;
1129   ptbins[2]=1.;
1130   ptbins[3]=2.;
1131   ptbins[4]=3.;
1132   ptbins[5]=4.;
1133   ptbins[6]=5.;
1134   ptbins[7]=6.;
1135   ptbins[8]=8.;
1136   ptbins[9]=12.;
1137   ptbins[10]=16.;
1138   ptbins[11]=20.;
1139   ptbins[12]=24.;
1140   ptbins[13]=ptmax;
1141
1142   SetGlobalIndex(nvars,nptbins);
1143   SetPtBins(nptbins+1,ptbins);
1144
1145   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73},/* pt<0.5*/
1146                                                   {0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73},/* 0.5<pt<1*/
1147                                                   {0.400,200.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.75},/* 1<pt<2 */
1148                                                   {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.8},/* 2<pt<3 */
1149                                                   {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 3<pt<4 */
1150                                                   {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 4<pt<5 */
1151                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 5<pt<6 */
1152                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85},/* 6<pt<8 */
1153                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85},/* 8<pt<12 */
1154                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.85},/* 12<pt<16 */
1155                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85},/* 16<pt<20 */
1156                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85},/* 20<pt<24 */
1157                                                   {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85}};/* pt>24 */
1158
1159
1160   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1161   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1162   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1163
1164   for (Int_t ibin=0;ibin<nptbins;ibin++){
1165     for (Int_t ivar = 0; ivar<nvars; ivar++){
1166       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1167     }
1168   }
1169
1170   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1171   SetUseSpecialCuts(kTRUE);
1172   SetRemoveDaughtersFromPrim(kTRUE);
1173
1174   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1175   delete [] cutsMatrixTransposeStand;
1176   cutsMatrixTransposeStand=NULL;
1177
1178   // PID SETTINGS
1179   AliAODPidHF* pidObj=new AliAODPidHF();
1180   //pidObj->SetName("pid4D0");
1181   Int_t mode=1;
1182   const Int_t nlims=2;
1183   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1184   Bool_t compat=kTRUE; //effective only for this mode
1185   Bool_t asym=kTRUE;
1186   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1187   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1188   pidObj->SetMatch(mode);
1189   pidObj->SetPLimit(plims,nlims);
1190   pidObj->SetSigma(sigmas);
1191   pidObj->SetCompat(compat);
1192   pidObj->SetTPC(kTRUE);
1193   pidObj->SetTOF(kTRUE);
1194
1195   SetPidHF(pidObj);
1196   SetUsePID(kTRUE);
1197   SetUseDefaultPID(kFALSE);
1198
1199
1200   PrintAll();
1201
1202   delete pidObj;
1203   pidObj=NULL;
1204
1205   return;
1206
1207 }
1208
1209
1210 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1211   //
1212   //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
1213   //... EVOLVING SOON
1214   //
1215
1216   SetName("D0toKpiCutsStandard");
1217   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1218
1219   // PILE UP REJECTION
1220   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1221
1222   // EVENT CUTS
1223   SetMinVtxContr(1);
1224
1225
1226   // TRACKS ON SINGLE TRACKS
1227   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1228   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1229   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1230   esdTrackCuts->SetRequireITSRefit(kTRUE);
1231   //  esdTrackCuts->SetMinNClustersITS(4);
1232   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1233   esdTrackCuts->SetMinDCAToVertexXY(0.);
1234   esdTrackCuts->SetEtaRange(-0.8,0.8);
1235   esdTrackCuts->SetPtRange(0.3,1.e10);
1236
1237   AddTrackCuts(esdTrackCuts);
1238
1239   const Int_t nptbins =13;
1240   const Double_t ptmax = 9999.;
1241   const Int_t nvars=9;
1242   Float_t ptbins[nptbins+1];
1243   ptbins[0]=0.;
1244   ptbins[1]=0.5;
1245   ptbins[2]=1.;
1246   ptbins[3]=2.;
1247   ptbins[4]=3.;
1248   ptbins[5]=4.;
1249   ptbins[6]=5.;
1250   ptbins[7]=6.;
1251   ptbins[8]=8.;
1252   ptbins[9]=12.;
1253   ptbins[10]=16.;
1254   ptbins[11]=20.;
1255   ptbins[12]=24.;
1256   ptbins[13]=ptmax;
1257
1258   SetGlobalIndex(nvars,nptbins);
1259   SetPtBins(nptbins+1,ptbins);
1260
1261   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.8},/* pt<0.5*/
1262                                                   {0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.8},/* 0.5<pt<1*/
1263                                                   {0.400,250.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-32000.*1E-8,0.8},/* 1<pt<2 */
1264                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-26000.*1E-8,0.94},/* 2<pt<3 */
1265                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1500.*1E-8,0.88},/* 3<pt<4 */
1266                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-1500.*1E-8,0.88},/* 4<pt<5 */
1267                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 5<pt<6 */
1268                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 6<pt<8 */
1269                                                   {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-10000.*1E-8,0.90},/* 8<pt<12 */
1270                                                   {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.90},/* 12<pt<16 */
1271                                                   {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85},/* 16<pt<20 */
1272                                                   {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.85},/* 20<pt<24 */
1273                                                   {0.400,350.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.82}};/* pt>24 */
1274
1275
1276   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1277   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1278   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1279
1280   for (Int_t ibin=0;ibin<nptbins;ibin++){
1281     for (Int_t ivar = 0; ivar<nvars; ivar++){
1282       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1283     }
1284   }
1285
1286   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1287   SetUseSpecialCuts(kTRUE);
1288   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1289   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1290   delete [] cutsMatrixTransposeStand;
1291   cutsMatrixTransposeStand=NULL;
1292
1293   // PID SETTINGS
1294   AliAODPidHF* pidObj=new AliAODPidHF();
1295   //pidObj->SetName("pid4D0");
1296   Int_t mode=1;
1297   const Int_t nlims=2;
1298   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1299   Bool_t compat=kTRUE; //effective only for this mode
1300   Bool_t asym=kTRUE;
1301   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1302   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1303   pidObj->SetMatch(mode);
1304   pidObj->SetPLimit(plims,nlims);
1305   pidObj->SetSigma(sigmas);
1306   pidObj->SetCompat(compat);
1307   pidObj->SetTPC(kTRUE);
1308   pidObj->SetTOF(kTRUE);
1309
1310   SetPidHF(pidObj);
1311   SetUsePID(kTRUE);
1312   SetUseDefaultPID(kTRUE);// TEMPORARY: PROTON EXCLUSION SET ONLY IN DEFAULT PID
1313
1314
1315   PrintAll();
1316
1317
1318   delete pidObj;
1319   pidObj=NULL;
1320
1321   return;
1322
1323 }