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