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