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