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