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   for(Int_t daught=0;daught<2;daught++){
 756     //Loop con prongs
 757     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
 758     if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0;
 759
 760     if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
 761       checkPIDInfo[daught]=kFALSE;
 762       continue;
 763     }
 764     Double_t pProng=aodtrack->P();
 765
 766     // identify kaon
 767     if(pProng<fmaxPtrackForPID){
 768       combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 769     }
 770     // identify pion
 771     if(pProng<fmaxPtrackForPID){
 772       if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
 773         combinedPID[daught][1]=0;
 774       }else{
 775         fPidHF->SetTOF(kFALSE);
 776         combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
 777         fPidHF->SetTOF(kTRUE);
 778         fPidHF->SetCompat(kTRUE);
 779       }
 780     }
 781
 782     if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
 783       isD0D0barPID[0]=0;
 784       isD0D0barPID[1]=0;
 785     }
 786     else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
 787       if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
 788       else isD0D0barPID[0]=0;// if K+ D0 excluded
 789     }
 790     /*    else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
 791           isD0D0barPID[0]=0;
 792           isD0D0barPID[1]=0;
 793           }
 794     */
 795     else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
 796       if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
 797       else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
 798     }
 799     else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
 800       if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
 801       else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
 802     }
 803
 804     if(fLowPt && d->Pt()<fPtLowPID){
 805      Double_t sigmaTPC[3]={3.,2.,0.};
 806      fPidHF->SetSigmaForTPC(sigmaTPC);
 807     // identify kaon
 808     combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 809
 810     if(pProng<0.6){
 811      fPidHF->SetCompat(kFALSE);
 812      combinedPID[daught][0]=fPidHF->MakeRawPid(aodtrack,3);
 813      fPidHF->SetCompat(kTRUE);
 814     }
 815
 816     if(!(fPidHF->CheckStatus(aodtrack,"TPC"))) {
 817      combinedPID[daught][1]=0;
 818     }else{
 819       fPidHF->SetTOF(kFALSE);
 820       Double_t sigmaTPCpi[3]={3.,3.,0.};
 821       fPidHF->SetSigmaForTPC(sigmaTPCpi);
 822       combinedPID[daught][1]=fPidHF->MakeRawPid(aodtrack,2);
 823       fPidHF->SetTOF(kTRUE);
 824        if(pProng<0.8){
 825         Bool_t isTPCpion=fPidHF->IsPionRaw(aodtrack,"TPC");
 826         if(isTPCpion){
 827          combinedPID[daught][1]=1;
 828         }else{
 829          combinedPID[daught][1]=-1;
 830         }
 831       }
 832     }
 833
 834    }
 835    fPidHF->SetSigmaForTPC(sigma_tmp);
 836   }// END OF LOOP ON DAUGHTERS
 837
 838    if(!checkPIDInfo[0] && !checkPIDInfo[1]) {
 839     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 840     return 0;
 841    }
 842
 843
 844   // FURTHER PID REQUEST (both daughter info is needed)
 845   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
 846     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
 847     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 848     return 0;
 849   }
 850
 851   if(fLowPt && d->Pt()<fPtLowPID){
 852     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
 853       fWhyRejection=32;// reject cases where the Kaon is not identified
 854       fPidHF->SetSigmaForTPC(sigma_tmp);
 855       return 0;
 856     }
 857   }
 858     if(fLowPt) fPidHF->SetSigmaForTPC(sigma_tmp);
 859
 860   //  cout<<"Why? "<<fWhyRejection<<endl;
 861   return isD0D0barPID[0]+isD0D0barPID[1];
 862 }
 863 //---------------------------------------------------------------------------
 864 Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
 865 {
 866   // ############################################################
 867   //
 868   // Apply PID selection
 869   //
 870   //
 871   // temporary selection: PID AS USED FOR D0 by Andrea Rossi (up to 28/06/2010)
 872   //
 873   // d must be a AliAODRecoDecayHF2Prong object
 874   // returns 0 if both D0 and D0bar are rejectecd
 875   //         1 if D0 is accepted while D0bar is rejected
 876   //         2 if D0bar is accepted while D0 is rejected
 877   //         3 if both are accepted
 878   // fWhyRejection variable (not returned for the moment, print it if needed)
 879   //               keeps some information on why a candidate has been
 880   //               rejected according to the following (unfriendly?) scheme
 881   //             if more rejection cases are considered interesting, just add numbers
 882   //
 883   //      TO BE CONSIDERED WITH A GRAIN OF SALT (the order in which cut are applied is relevant)
 884   //              from 20 to 30: "detector" selection (PID acceptance)
 885   //                                                  26: TPC refit
 886   //                                                  27: ITS refit
 887   //                                                  28: no (TOF||TPC) pid information (no kTOFpid,kTOFout,kTIME,kTPCpid,...)
 888   //
 889   //              from 30 to 40: PID selection
 890   //                                                  31: no Kaon compatible tracks found between daughters
 891   //                                                  32: no Kaon identified tracks found (strong sel. at low momenta)
 892   //                                                  33: both mass hypotheses are rejected
 893   //
 894   // ############################################################
 895
 896   if(!fUsePID) return 3;
 897   fWhyRejection=0;
 898   Int_t isD0D0barPID[2]={1,2};
 899   Double_t nsigmaTPCpi=-1., nsigmaTPCK=-1.; //used for TPC pid
 900   Double_t tofSig,times[5];// used fot TOF pid
 901   Int_t hasPID[2]={2,2};// flag to count how many detectors give PID info for the daughters
 902   Int_t isKaonPionTOF[2][2],isKaonPionTPC[2][2];
 903   Int_t combinedPID[2][2];// CONVENTION: [daught][isK,IsPi]; [0][0]=(prong 1, isK)=value [0][1]=(prong 1, isPi)=value;
 904   //                                                                                                 same for prong 2
 905   //                                               values convention -1 = discarded
 906   //                                                                  0 = not identified (but compatible) || No PID (->hasPID flag)
 907   //                                                                  1 = identified
 908   // PID search:   pion (TPC) or not K (TOF), Kaon hypothesis for both
 909   // Initial hypothesis: unknwon (but compatible)
 910   isKaonPionTOF[0][0]=0;
 911   isKaonPionTOF[0][1]=0;
 912   isKaonPionTOF[1][0]=0;
 913   isKaonPionTOF[1][1]=0;
 914
 915   isKaonPionTPC[0][0]=0;
 916   isKaonPionTPC[0][1]=0;
 917   isKaonPionTPC[1][0]=0;
 918   isKaonPionTPC[1][1]=0;
 919
 920   combinedPID[0][0]=0;
 921   combinedPID[0][1]=0;
 922   combinedPID[1][0]=0;
 923   combinedPID[1][1]=0;
 924
 925
 926
 927   for(Int_t daught=0;daught<2;daught++){
 928     //Loop con prongs
 929
 930     // ########### Step 0- CHECKING minimal PID "ACCEPTANCE" ####################
 931
 932     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
 933
 934     if(!(aodtrack->GetStatus()&AliESDtrack::kTPCrefit)){
 935       fWhyRejection=26;
 936       return 0;
 937     }
 938     if(!(aodtrack->GetStatus()&AliESDtrack::kITSrefit)){
 939       fWhyRejection=27;
 940       return 0;
 941     }
 942
 943     AliAODPid *pid=aodtrack->GetDetPid();
 944     if(!pid) {
 945       //delete esdtrack;
 946       hasPID[daught]--;
 947       continue;
 948     }
 949
 950     // ########### Step 1- Check of TPC and TOF response ####################
 951
 952     Double_t ptrack=aodtrack->P();
 953     //#################### TPC PID #######################
 954      if (!(aodtrack->GetStatus()&AliESDtrack::kTPCpid )){
 955        // NO TPC PID INFO FOR THIS TRACK
 956        hasPID[daught]--;
 957      }
 958      else {
 959        static AliTPCPIDResponse theTPCpid;
 960        AliAODPid *pidObj = aodtrack->GetDetPid();
 961        Double_t ptProng=pidObj->GetTPCmomentum();
 962        nsigmaTPCpi = theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kPion);
 963        nsigmaTPCK =  theTPCpid.GetNumberOfSigmas(ptProng,(Float_t)pid->GetTPCsignal(),(Int_t)aodtrack->GetTPCClusterMap().CountBits(),AliPID::kKaon);
 964        //if(ptrack<0.6){
 965        if(ptProng<0.6){
 966          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
 967          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 968          if(TMath::Abs(nsigmaTPCpi)<2.)isKaonPionTPC[daught][1]=1;
 969          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 970        }
 971        //else if(ptrack<.8){
 972        else if(ptProng<.8){
 973          if(TMath::Abs(nsigmaTPCK)<1.)isKaonPionTPC[daught][0]=1;
 974          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 975          if(TMath::Abs(nsigmaTPCpi)<1.)isKaonPionTPC[daught][1]=1;
 976          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 977        }
 978        else {
 979          //     if(nsigmaTPCK>-2.&&nsigmaTPCK<1.)isKaonPionTPC[daught][0]=1;
 980          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
 981          //if(nsigmaTPCpi>-1.&&nsigmaTPCpi<2.)isKaonPionTPC[daught][1]=1;
 982          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
 983        }
 984      }
 985
 986
 987     // ##### TOF PID: do not ask nothing for pion/protons ############
 988      if(!((aodtrack->GetStatus()&AliESDtrack::kTOFpid)&&(aodtrack->GetStatus()&AliESDtrack::kTOFout)&&(aodtrack->GetStatus()&AliESDtrack::kTIME))){
 989        // NO TOF PID INFO FOR THIS TRACK
 990        hasPID[daught]--;
 991      }
 992      else{
 993        tofSig=pid->GetTOFsignal();
 994        pid->GetIntegratedTimes(times);
 995        if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
 996        if(TMath::Abs(tofSig-times[3])>3.*160.){
 997          isKaonPionTOF[daught][0]=-1;
 998        }
 999        else {
1000          if(ptrack<1.5){
1001            isKaonPionTOF[daught][0]=1;
1002          }
1003        }
1004      }
1005
1006      //######### Step 2: COMBINE TOF and TPC PID ###############
1007      // we apply the following convention: if TPC and TOF disagree (discarded Vs identified) -> unknown
1008      combinedPID[daught][0]=isKaonPionTOF[daught][0]+isKaonPionTPC[daught][0];
1009      combinedPID[daught][1]=isKaonPionTOF[daught][1]+isKaonPionTPC[daught][1];
1010
1011
1012      //######### Step 3:   USE PID INFO
1013
1014      if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){// if not a K- and not a pi- both D0 and D0bar excluded
1015        isD0D0barPID[0]=0;
1016        isD0D0barPID[1]=0;
1017      }
1018      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
1019        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
1020        else isD0D0barPID[0]=0;// if K+ D0 excluded
1021      }
1022      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
1023        isD0D0barPID[0]=0;
1024        isD0D0barPID[1]=0;
1025      }
1026      else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){
1027        if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
1028        else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
1029      }
1030      else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
1031        if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
1032       else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
1033      }
1034
1035      // ##########  ALSO DIFFERENT TPC PID REQUEST FOR LOW pt D0: request of K identification      ###############################
1036      // ########## more tolerant criteria for single particle ID-> more selective criteria for D0   ##############################
1037      // ###############                     NOT OPTIMIZED YET                                  ###################################
1038      if(d->Pt()<2.){
1039        isKaonPionTPC[daught][0]=0;
1040        isKaonPionTPC[daught][1]=0;
1041        AliAODPid *pidObj = aodtrack->GetDetPid();
1042        Double_t ptProng=pidObj->GetTPCmomentum();
1043        //if(ptrack<0.6){
1044        if(ptProng<0.6){
1045          if(TMath::Abs(nsigmaTPCK)<3.)isKaonPionTPC[daught][0]=1;
1046          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1047          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1048          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1049      }
1050        //else if(ptrack<.8){
1051        else if(ptProng<.8){
1052          if(TMath::Abs(nsigmaTPCK)<2.)isKaonPionTPC[daught][0]=1;
1053          else if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1054          if(TMath::Abs(nsigmaTPCpi)<3.)isKaonPionTPC[daught][1]=1;
1055          else if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1056        }
1057        else {
1058          if(TMath::Abs(nsigmaTPCK)>3.)isKaonPionTPC[daught][0]=-1;
1059          if(TMath::Abs(nsigmaTPCpi)>3.)isKaonPionTPC[daught][1]=-1;
1060        }
1061      }
1062
1063   }// END OF LOOP ON DAUGHTERS
1064
1065   // FURTHER PID REQUEST (both daughter info is needed)
1066   if(combinedPID[0][0]<=-1&&combinedPID[1][0]<=-1){
1067     fWhyRejection=31;// reject cases in which no kaon-compatible tracks are found
1068     return 0;
1069   }
1070   else if(hasPID[0]==0&&hasPID[1]==0){
1071     fWhyRejection=28;// reject cases in which no PID info is available
1072     return 0;
1073   }
1074   if(d->Pt()<2.){
1075     // request of K identification at low D0 pt
1076     combinedPID[0][0]=0;
1077     combinedPID[0][1]=0;
1078     combinedPID[1][0]=0;
1079     combinedPID[1][1]=0;
1080
1081     combinedPID[0][0]=isKaonPionTOF[0][0]+isKaonPionTPC[0][0];
1082     combinedPID[0][1]=isKaonPionTOF[0][1]+isKaonPionTPC[0][1];
1083     combinedPID[1][0]=isKaonPionTOF[1][0]+isKaonPionTPC[1][0];
1084     combinedPID[1][1]=isKaonPionTOF[1][1]+isKaonPionTPC[1][1];
1085
1086     if(combinedPID[0][0]<=0&&combinedPID[1][0]<=0){
1087       fWhyRejection=32;// reject cases where the Kaon is not identified
1088       return 0;
1089     }
1090   }
1091
1092   //  cout<<"Why? "<<fWhyRejection<<endl;
1093   return isD0D0barPID[0]+isD0D0barPID[1];
1094 }
1095
1096
1097
1098 //---------------------------------------------------------------------------
1099 Int_t AliRDHFCutsD0toKpi::CombineSelectionLevels(Int_t selectionvalTrack,
1100                                                  Int_t selectionvalCand,
1101                                                  Int_t selectionvalPID) const
1102 {
1103   //
1104   // This method combines the tracks, PID and cuts selection results
1105   //
1106   if(selectionvalTrack==0) return 0;
1107
1108   Int_t returnvalue;
1109
1110   switch(selectionvalPID) {
1111   case 0:
1112     returnvalue=0;
1113     break;
1114   case 1:
1115     returnvalue=((selectionvalCand==1 || selectionvalCand==3) ? 1 : 0);
1116     break;
1117   case 2:
1118     returnvalue=((selectionvalCand==2 || selectionvalCand==3) ? 2 : 0);
1119     break;
1120   case 3:
1121     returnvalue=selectionvalCand;
1122     break;
1123   default:
1124     returnvalue=0;
1125     break;
1126   }
1127
1128   return returnvalue;
1129 }
1130 //----------------------------------------------------------------------------
1131 Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
1132 {
1133   //
1134   // Note: this method is temporary
1135   // Additional cuts on decay lenght and lower cut for d0 norm are applied using vertex without candidate's daughters
1136   //
1137
1138   //apply cuts
1139
1140   Float_t normDecLengthCut=1.,decLengthCut=TMath::Min(d->P()*0.0066+0.01,0.06/*cm*/), normd0Cut=0.5;
1141   // "decay length" expo law with tau' = beta*gamma*ctau= p/m*ctau =p*0.0123/1.864~p*0.0066
1142   // decay lenght > ctau' implies to retain (1-1/e) (for signal without considering detector resolution),
1143
1144   Int_t returnvalue=3; //cut passed
1145   for(Int_t i=0;i<2/*prongs*/;i++){
1146     if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
1147   }
1148   if(d->DecayLength2()<decLengthCut*decLengthCut)  return 0; //decLengthCut not passed
1149   if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut)  return 0; //decLengthCut not passed
1150
1151   return returnvalue;
1152 }
1153
1154 //----------------------------------------------
1155 void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
1156 {
1157   //switch on candidate selection via AliKFparticle
1158   if(!useKF) return;
1159   if(useKF){
1160     fUseKF=useKF;
1161     Int_t nvarsKF=11;
1162     SetNVars(nvarsKF);
1163     TString varNamesKF[11]={"inv. mass [GeV]",
1164                             "dca [cm]",
1165                             "cosThetaStar",
1166                             "pTK [GeV/c]",
1167                             "pTPi [GeV/c]",
1168                             "d0K [cm]",
1169                             "d0Pi [cm]",
1170                             "d0d0 [cm^2]",
1171                             "cosThetaPoint"
1172                             "DecayLength[cm]",
1173                             "RedChi2"};
1174     Bool_t isUpperCutKF[11]={kTRUE,
1175                              kTRUE,
1176                              kTRUE,
1177                              kFALSE,
1178                              kFALSE,
1179                              kTRUE,
1180                              kTRUE,
1181                              kTRUE,
1182                              kFALSE,
1183                              kFALSE,
1184                              kTRUE};
1185     SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
1186     SetGlobalIndex();
1187     Bool_t forOpt[11]={kFALSE,
1188                     kTRUE,
1189                     kTRUE,
1190                     kFALSE,
1191                     kFALSE,
1192                     kFALSE,
1193                     kFALSE,
1194                     kTRUE,
1195                     kTRUE,
1196                     kFALSE,
1197                     kFALSE};
1198     SetVarsForOpt(4,forOpt);
1199   }
1200   return;
1201 }
1202
1203
1204 //---------------------------------------------------------------------------
1205 void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
1206   //
1207   //STANDARD CUTS USED FOR 2010 pp analysis
1208   //dca cut will be enlarged soon to 400 micron
1209   //
1210
1211   SetName("D0toKpiCutsStandard");
1212   SetTitle("Standard Cuts for D0 analysis");
1213
1214   // PILE UP REJECTION
1215   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1216
1217   // EVENT CUTS
1218   SetMinVtxContr(1);
1219  // MAX Z-VERTEX CUT
1220   SetMaxVtxZ(10.);
1221
1222   // TRACKS ON SINGLE TRACKS
1223   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1224   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1225   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1226   esdTrackCuts->SetRequireITSRefit(kTRUE);
1227   //  esdTrackCuts->SetMinNClustersITS(4);
1228   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1229   esdTrackCuts->SetMinDCAToVertexXY(0.);
1230   esdTrackCuts->SetEtaRange(-0.8,0.8);
1231   esdTrackCuts->SetPtRange(0.3,1.e10);
1232
1233   AddTrackCuts(esdTrackCuts);
1234
1235
1236   const Int_t nptbins =14;
1237   const Double_t ptmax = 9999.;
1238   const Int_t nvars=11;
1239   Float_t ptbins[nptbins+1];
1240   ptbins[0]=0.;
1241   ptbins[1]=0.5;
1242   ptbins[2]=1.;
1243   ptbins[3]=2.;
1244   ptbins[4]=3.;
1245   ptbins[5]=4.;
1246   ptbins[6]=5.;
1247   ptbins[7]=6.;
1248   ptbins[8]=7.;
1249   ptbins[9]=8.;
1250   ptbins[10]=12.;
1251   ptbins[11]=16.;
1252   ptbins[12]=20.;
1253   ptbins[13]=24.;
1254   ptbins[14]=ptmax;
1255
1256   SetGlobalIndex(nvars,nptbins);
1257   SetPtBins(nptbins+1,ptbins);
1258
1259   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*/
1260                                                   {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*/
1261                                                   {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 */
1262                                                   {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 */
1263                                                   {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 */
1264                                                   {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 */
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.},/* 5<pt<6 */
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.},/* 6<pt<7 */
1267                                                   {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 */
1268                                                   {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 */
1269                                                   {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 */
1270                                                   {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 */
1271                                                   {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 */
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.}};/* pt>24 */
1273
1274
1275   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1276   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1277   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1278
1279   for (Int_t ibin=0;ibin<nptbins;ibin++){
1280     for (Int_t ivar = 0; ivar<nvars; ivar++){
1281       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1282     }
1283   }
1284
1285   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1286   SetUseSpecialCuts(kTRUE);
1287   SetRemoveDaughtersFromPrim(kTRUE);
1288
1289   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1290   delete [] cutsMatrixTransposeStand;
1291   cutsMatrixTransposeStand=NULL;
1292
1293   // PID SETTINGS
1294   AliAODPidHF* pidObj=new AliAODPidHF();
1295   //pidObj->SetName("pid4D0");
1296   Int_t mode=1;
1297   const Int_t nlims=2;
1298   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1299   Bool_t compat=kTRUE; //effective only for this mode
1300   Bool_t asym=kTRUE;
1301   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1302   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1303   pidObj->SetMatch(mode);
1304   pidObj->SetPLimit(plims,nlims);
1305   pidObj->SetSigma(sigmas);
1306   pidObj->SetCompat(compat);
1307   pidObj->SetTPC(kTRUE);
1308   pidObj->SetTOF(kTRUE);
1309   pidObj->SetPCompatTOF(1.5);
1310   pidObj->SetSigmaForTPCCompat(3.);
1311   pidObj->SetSigmaForTOFCompat(3.);
1312   pidObj->SetOldPid(kTRUE);
1313
1314   SetPidHF(pidObj);
1315   SetUsePID(kTRUE);
1316   SetUseDefaultPID(kFALSE);
1317
1318   SetLowPt(kFALSE);
1319
1320   PrintAll();
1321
1322   delete pidObj;
1323   pidObj=NULL;
1324
1325   return;
1326
1327 }
1328
1329
1330 //---------------------------------------------------------------------------
1331 void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
1332   //
1333   // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1334   //
1335
1336   SetName("D0toKpiCutsStandard");
1337   SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
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   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1355   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1356   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1357
1358   AddTrackCuts(esdTrackCuts);
1359
1360
1361   const Int_t nvars=11;
1362   const Int_t nptbins=13;
1363   Float_t ptbins[nptbins+1];
1364   ptbins[0]=0.;
1365   ptbins[1]=0.5;
1366   ptbins[2]=1.;
1367   ptbins[3]=2.;
1368   ptbins[4]=3.;
1369   ptbins[5]=4.;
1370   ptbins[6]=5.;
1371   ptbins[7]=6.;
1372   ptbins[8]=8.;
1373   ptbins[9]=12.;
1374   ptbins[10]=16.;
1375   ptbins[11]=20.;
1376   ptbins[12]=24.;
1377   ptbins[13]=9999.;
1378
1379   SetPtBins(nptbins+1,ptbins);
1380
1381
1382   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*/
1383                                                   {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*/
1384                                                   {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 */
1385                                                   {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 */
1386                                                   {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 */
1387                                                   {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 */
1388                                                   {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 */
1389                                                   {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1390                                                   {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 */
1391                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1392                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1393                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1394                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1395
1396   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1397   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1398   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1399   for (Int_t ibin=0;ibin<nptbins;ibin++){
1400     for (Int_t ivar = 0; ivar<nvars; ivar++){
1401       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1402     }
1403   }
1404   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1405   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1406   delete [] cutsMatrixTransposeStand;
1407
1408
1409   //pid settings
1410   AliAODPidHF* pidObj=new AliAODPidHF();
1411   Int_t mode=1;
1412   const Int_t nlims=2;
1413   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1414   Bool_t compat=kTRUE; //effective only for this mode
1415   Bool_t asym=kTRUE;
1416   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1417   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1418   pidObj->SetMatch(mode);
1419   pidObj->SetPLimit(plims,nlims);
1420   pidObj->SetSigma(sigmas);
1421   pidObj->SetCompat(compat);
1422   pidObj->SetTPC(kTRUE);
1423   pidObj->SetTOF(kTRUE);
1424   pidObj->SetOldPid(kTRUE);
1425   SetPidHF(pidObj);
1426
1427   SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1428
1429   SetUsePID(kTRUE);
1430
1431   //activate pileup rejection (for pp)
1432   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1433
1434   //Do recalculate the vertex
1435   SetRemoveDaughtersFromPrim(kTRUE);
1436
1437   // Cut on the zvtx
1438   SetMaxVtxZ(10.);
1439
1440   // Use the kFirst selection for tracks with candidate D with pt<2
1441   SetSelectCandTrackSPDFirst(kTRUE,2.);
1442
1443   // Use special cuts for D candidates with pt<2
1444   SetUseSpecialCuts(kTRUE);
1445   SetMaximumPtSpecialCuts(2.);
1446
1447   PrintAll();
1448
1449   delete pidObj;
1450   pidObj=NULL;
1451
1452   return;
1453 }
1454
1455
1456 //---------------------------------------------------------------------------
1457 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1458   //
1459   // Final CUTS USED FOR 2010 PbPb analysis of central events
1460   // REMEMBER TO EVENTUALLY SWITCH ON
1461   //          SetUseAOD049(kTRUE);
1462   //
1463
1464   SetName("D0toKpiCutsStandard");
1465   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1466
1467
1468   // CENTRALITY SELECTION
1469   SetMinCentrality(0.);
1470   SetMaxCentrality(80.);
1471   SetUseCentrality(AliRDHFCuts::kCentV0M);
1472
1473
1474
1475   // EVENT CUTS
1476   SetMinVtxContr(1);
1477   // MAX Z-VERTEX CUT
1478   SetMaxVtxZ(10.);
1479   // PILE UP
1480   SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1481   // PILE UP REJECTION
1482   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1483
1484   // TRACKS ON SINGLE TRACKS
1485   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1486   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1487   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1488   esdTrackCuts->SetRequireITSRefit(kTRUE);
1489   //  esdTrackCuts->SetMinNClustersITS(4);
1490   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1491   esdTrackCuts->SetMinDCAToVertexXY(0.);
1492   esdTrackCuts->SetEtaRange(-0.8,0.8);
1493   esdTrackCuts->SetPtRange(0.7,1.e10);
1494
1495   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1496   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1497   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1498
1499
1500   AddTrackCuts(esdTrackCuts);
1501   // SPD k FIRST for Pt<3 GeV/c
1502   SetSelectCandTrackSPDFirst(kTRUE, 3);
1503
1504   // CANDIDATE CUTS
1505   const Int_t nptbins =13;
1506   const Double_t ptmax = 9999.;
1507   const Int_t nvars=11;
1508   Float_t ptbins[nptbins+1];
1509   ptbins[0]=0.;
1510   ptbins[1]=0.5;
1511   ptbins[2]=1.;
1512   ptbins[3]=2.;
1513   ptbins[4]=3.;
1514   ptbins[5]=4.;
1515   ptbins[6]=5.;
1516   ptbins[7]=6.;
1517   ptbins[8]=8.;
1518   ptbins[9]=12.;
1519   ptbins[10]=16.;
1520   ptbins[11]=20.;
1521   ptbins[12]=24.;
1522   ptbins[13]=ptmax;
1523
1524   SetGlobalIndex(nvars,nptbins);
1525   SetPtBins(nptbins+1,ptbins);
1526   SetMinPtCandidate(2.);
1527
1528   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*/
1529                                                   {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*/
1530                                                   {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 */
1531                                                   {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 */
1532                                                   {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 */
1533                                                   {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 */
1534                                                   {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 */
1535                                                   {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 */
1536                                                   {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 */
1537                                                   {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 */
1538                                                   {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 */
1539                                                   {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 */
1540                                                   {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 */
1541
1542
1543   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1544   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1545   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1546
1547   for (Int_t ibin=0;ibin<nptbins;ibin++){
1548     for (Int_t ivar = 0; ivar<nvars; ivar++){
1549       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1550     }
1551   }
1552
1553   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1554   SetUseSpecialCuts(kTRUE);
1555   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1556   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1557   delete [] cutsMatrixTransposeStand;
1558   cutsMatrixTransposeStand=NULL;
1559
1560   // PID SETTINGS
1561   AliAODPidHF* pidObj=new AliAODPidHF();
1562   //pidObj->SetName("pid4D0");
1563   Int_t mode=1;
1564   const Int_t nlims=2;
1565   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1566   Bool_t compat=kTRUE; //effective only for this mode
1567   Bool_t asym=kTRUE;
1568   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1569   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1570   pidObj->SetMatch(mode);
1571   pidObj->SetPLimit(plims,nlims);
1572   pidObj->SetSigma(sigmas);
1573   pidObj->SetCompat(compat);
1574   pidObj->SetTPC(kTRUE);
1575   pidObj->SetTOF(kTRUE);
1576   pidObj->SetPCompatTOF(2.);
1577   pidObj->SetSigmaForTPCCompat(3.);
1578   pidObj->SetSigmaForTOFCompat(3.);
1579   pidObj->SetOldPid(kTRUE);
1580
1581
1582   SetPidHF(pidObj);
1583   SetUsePID(kTRUE);
1584   SetUseDefaultPID(kFALSE);
1585
1586
1587   PrintAll();
1588
1589
1590   delete pidObj;
1591   pidObj=NULL;
1592
1593   return;
1594
1595 }
1596
1597 //---------------------------------------------------------------------------
1598 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
1599   // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1600
1601
1602   SetName("D0toKpiCutsStandard");
1603   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1604
1605
1606   // CENTRALITY SELECTION
1607   SetMinCentrality(40.);
1608   SetMaxCentrality(80.);
1609   SetUseCentrality(AliRDHFCuts::kCentV0M);
1610
1611   // EVENT CUTS
1612   SetMinVtxContr(1);
1613   // MAX Z-VERTEX CUT
1614   SetMaxVtxZ(10.);
1615   // PILE UP
1616   SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1617   // PILE UP REJECTION
1618   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1619
1620   // TRACKS ON SINGLE TRACKS
1621   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1622   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1623   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1624   esdTrackCuts->SetRequireITSRefit(kTRUE);
1625   //  esdTrackCuts->SetMinNClustersITS(4);
1626   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1627   esdTrackCuts->SetMinDCAToVertexXY(0.);
1628   esdTrackCuts->SetEtaRange(-0.8,0.8);
1629   esdTrackCuts->SetPtRange(0.5,1.e10);
1630
1631   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1632   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1633   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1634
1635
1636   AddTrackCuts(esdTrackCuts);
1637   // SPD k FIRST for Pt<3 GeV/c
1638   SetSelectCandTrackSPDFirst(kTRUE, 3);
1639
1640   // CANDIDATE CUTS
1641   const Int_t nptbins =13;
1642   const Double_t ptmax = 9999.;
1643   const Int_t nvars=11;
1644   Float_t ptbins[nptbins+1];
1645   ptbins[0]=0.;
1646   ptbins[1]=0.5;
1647   ptbins[2]=1.;
1648   ptbins[3]=2.;
1649   ptbins[4]=3.;
1650   ptbins[5]=4.;
1651   ptbins[6]=5.;
1652   ptbins[7]=6.;
1653   ptbins[8]=8.;
1654   ptbins[9]=12.;
1655   ptbins[10]=16.;
1656   ptbins[11]=20.;
1657   ptbins[12]=24.;
1658   ptbins[13]=ptmax;
1659
1660   SetGlobalIndex(nvars,nptbins);
1661   SetPtBins(nptbins+1,ptbins);
1662   SetMinPtCandidate(0.);
1663
1664   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*/
1665                                                   {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*/
1666                                                   {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 */
1667                                                   {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 */
1668                                                   {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 */
1669                                                   {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 */
1670                                                   {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 */
1671                                                   {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 */
1672                                                   {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 */
1673                                                   {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 */
1674                                                   {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 */
1675                                                   {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 */
1676                                                   {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 */
1677
1678
1679   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1680   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1681   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1682
1683   for (Int_t ibin=0;ibin<nptbins;ibin++){
1684     for (Int_t ivar = 0; ivar<nvars; ivar++){
1685       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1686     }
1687   }
1688
1689   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1690   SetUseSpecialCuts(kTRUE);
1691   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1692   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1693   delete [] cutsMatrixTransposeStand;
1694   cutsMatrixTransposeStand=NULL;
1695
1696   // PID SETTINGS
1697   AliAODPidHF* pidObj=new AliAODPidHF();
1698   //pidObj->SetName("pid4D0");
1699   Int_t mode=1;
1700   const Int_t nlims=2;
1701   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1702   Bool_t compat=kTRUE; //effective only for this mode
1703   Bool_t asym=kTRUE;
1704   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1705   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1706   pidObj->SetMatch(mode);
1707   pidObj->SetPLimit(plims,nlims);
1708   pidObj->SetSigma(sigmas);
1709   pidObj->SetCompat(compat);
1710   pidObj->SetTPC(kTRUE);
1711   pidObj->SetTOF(kTRUE);
1712   pidObj->SetPCompatTOF(2.);
1713   pidObj->SetSigmaForTPCCompat(3.);
1714   pidObj->SetSigmaForTOFCompat(3.);
1715   pidObj->SetOldPid(kTRUE);
1716
1717   SetPidHF(pidObj);
1718   SetUsePID(kTRUE);
1719   SetUseDefaultPID(kFALSE);
1720
1721   SetLowPt(kTRUE,2.);
1722
1723   PrintAll();
1724
1725
1726   delete pidObj;
1727   pidObj=NULL;
1728
1729   return;
1730
1731 }
1732
1733
1734 //---------------------------------------------------------------------------
1735 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
1736
1737   // Default 2010 PbPb cut object
1738   SetStandardCutsPbPb2010();
1739   AliAODPidHF *pidobj=GetPidHF();
1740
1741   pidobj->SetOldPid(kFALSE);
1742
1743   //
1744   // Enable all 2011 PbPb run triggers
1745   //
1746   SetTriggerClass("");
1747   ResetMaskAndEnableMBTrigger();
1748   EnableCentralTrigger();
1749   EnableSemiCentralTrigger();
1750
1751 }
1752
1753 //---------------------------------------------------------------------------
1754 Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
1755 {
1756   // ############################################################
1757   //
1758   // Apply Bayesian PID selection
1759   // Implementation of Bayesian PID: Jeremy Wilkinson
1760   //
1761   // ############################################################
1762
1763   if (!fUsePID || !d) return 3;
1764
1765
1766   if (fBayesianStrategy == kBayesWeightNoFilter) {
1767      //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1768      CalculateBayesianWeights(d);
1769      return 3;
1770
1771 }
1772
1773
1774
1775   Int_t isD0 = 0;
1776   Int_t isD0bar = 0;
1777   Int_t returnvalue = 0;
1778
1779   Int_t isPosKaon = 0, isNegKaon = 0;
1780
1781   //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1782
1783
1784   Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
1785   AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1786
1787   if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
1788     return 0;  //reject if charges not opposite
1789   }
1790    Double_t momentumpositive=0., momentumnegative=0.;   //Used in "prob > prior" method
1791   for (Int_t daught = 0; daught < 2; daught++) {
1792     //Loop over prongs
1793
1794       if (aodtrack[daught]->Charge() == -1) {
1795          momentumnegative = aodtrack[daught]->P();
1796       }
1797       if (aodtrack[daught]->Charge() == +1) {
1798          momentumpositive = aodtrack[daught]->P();
1799       }
1800           if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1801              checkPIDInfo[daught] = kFALSE;
1802              continue;
1803           }
1804
1805    }
1806
1807    //Loop over daughters ends here
1808
1809    if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
1810       return 0;      //Reject if both daughters lack both TPC and TOF info
1811    }
1812
1813
1814    CalculateBayesianWeights(d);        //Calculates all Bayesian probabilities for both positive and negative tracks
1815     //      Double_t prob0[AliPID::kSPECIES];
1816     //      fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
1817    ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
1818    switch (fBayesianCondition) {
1819       ///A: Standard max. probability method (accept most likely species)
1820       case kMaxProb:
1821                         if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
1822                           isPosKaon = 1;  //flag [daught] as a kaon
1823                         }
1824
1825                         if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
1826                           isNegKaon = 1;  //flag [daught] as a kaon
1827                         }
1828                     break;
1829       ///B: Accept if probability greater than prior
1830       case kAbovePrior:
1831
1832                         if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1833                                                                         GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) {  //Retrieves relevant prior, gets value at momentum
1834                         isNegKaon = 1;
1835                         }
1836                         if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1837                                                                         GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) {  //Retrieves relevant prior, gets value at momentum
1838                         isPosKaon = 1;
1839                         }
1840
1841                 break;
1842
1843       ///C: Accept if probability greater than user-defined threshold
1844       case kThreshold:
1845          if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
1846             isNegKaon = 1;
1847          }
1848          if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
1849             isPosKaon = 1;
1850          }
1851
1852          break;
1853   }
1854
1855
1856      //Momentum-based selection (also applied to filtered weighted method)
1857
1858      if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
1859          if (isNegKaon && isPosKaon) { // If both are kaons, reject
1860             isD0 = 0;
1861             isD0bar = 0;
1862          } else if (isNegKaon) {       //If negative kaon present, D0
1863             isD0 = 1;
1864          } else if (isPosKaon) {       //If positive kaon present, D0bar
1865             isD0bar = 1;
1866          } else {                      //If neither ID'd as kaon, subject to extra tests
1867             isD0 = 1;
1868             isD0bar = 1;
1869             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
1870                if (aodtrack[0]->Charge() == -1) {
1871             isD0 = 0;  //Check charge and reject based on pion hypothesis
1872                }
1873                if (aodtrack[0]->Charge() == 1) {
1874             isD0bar = 0;
1875                }
1876             }
1877             if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
1878                if (aodtrack[1]->Charge() == -1) {
1879             isD0 = 0;
1880                }
1881                if (aodtrack[1]->Charge() == 1) {
1882             isD0bar = 0;
1883                }
1884             }
1885          }
1886
1887
1888
1889          if (isD0 && isD0bar) {
1890             returnvalue = 3;
1891          }
1892          if (isD0 && !isD0bar) {
1893             returnvalue = 1;
1894          }
1895          if (!isD0 && isD0bar) {
1896             returnvalue = 2;
1897          }
1898          if (!isD0 && !isD0bar) {
1899             returnvalue = 0;
1900          }
1901      }
1902
1903     //Simple Bayesian
1904     if (fBayesianStrategy == kBayesSimple) {
1905
1906          if (isPosKaon && isNegKaon)   {  //If both are ID'd as kaons, reject
1907                returnvalue = 0;
1908             } else if (isNegKaon)   {     //If negative kaon, D0
1909                returnvalue = 1;
1910             } else if (isPosKaon)   {     //If positive kaon, D0-bar
1911                returnvalue = 2;
1912             } else {
1913                returnvalue = 0;  //If neither kaon, reject
1914             }
1915     }
1916
1917   return returnvalue;
1918
1919
1920
1921 }
1922
1923
1924
1925 //---------------------------------------------------------------------------
1926 void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
1927
1928 {
1929   //Function to compute weights for Bayesian method
1930
1931   AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1932   if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
1933     return;  //Reject if charges do not oppose
1934   }
1935   for (Int_t daught = 0; daught < 2; daught++) {
1936     //Loop over prongs
1937
1938     if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1939       continue;
1940     }
1941
1942
1943     // identify kaon, define weights
1944     if (aodtrack[daught]->Charge() == +1) {
1945     fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
1946     }
1947
1948     if (aodtrack[daught]->Charge() == -1) {
1949     fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
1950     }
1951   }
1952 }
1953
1954
1955