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
1236
1237   const Int_t nptbins =14;
1238   const Double_t ptmax = 9999.;
1239   const Int_t nvars=11;
1240   Float_t ptbins[nptbins+1];
1241   ptbins[0]=0.;
1242   ptbins[1]=0.5;
1243   ptbins[2]=1.;
1244   ptbins[3]=2.;
1245   ptbins[4]=3.;
1246   ptbins[5]=4.;
1247   ptbins[6]=5.;
1248   ptbins[7]=6.;
1249   ptbins[8]=7.;
1250   ptbins[9]=8.;
1251   ptbins[10]=12.;
1252   ptbins[11]=16.;
1253   ptbins[12]=20.;
1254   ptbins[13]=24.;
1255   ptbins[14]=ptmax;
1256
1257   SetGlobalIndex(nvars,nptbins);
1258   SetPtBins(nptbins+1,ptbins);
1259
1260   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*/
1261                                                   {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*/
1262                                                   {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 */
1263                                                   {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 */
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.},/* 3<pt<4 */
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.},/* 4<pt<5 */
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.},/* 5<pt<6 */
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.},/* 6<pt<7 */
1268                                                   {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 */
1269                                                   {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 */
1270                                                   {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 */
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.},/* 16<pt<20 */
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.},/* 20<pt<24 */
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.}};/* pt>24 */
1274
1275
1276   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1277   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1278   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1279
1280   for (Int_t ibin=0;ibin<nptbins;ibin++){
1281     for (Int_t ivar = 0; ivar<nvars; ivar++){
1282       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1283     }
1284   }
1285
1286   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1287   SetUseSpecialCuts(kTRUE);
1288   SetRemoveDaughtersFromPrim(kTRUE);
1289
1290   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1291   delete [] cutsMatrixTransposeStand;
1292   cutsMatrixTransposeStand=NULL;
1293
1294   // PID SETTINGS
1295   AliAODPidHF* pidObj=new AliAODPidHF();
1296   //pidObj->SetName("pid4D0");
1297   Int_t mode=1;
1298   const Int_t nlims=2;
1299   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1300   Bool_t compat=kTRUE; //effective only for this mode
1301   Bool_t asym=kTRUE;
1302   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1303   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1304   pidObj->SetMatch(mode);
1305   pidObj->SetPLimit(plims,nlims);
1306   pidObj->SetSigma(sigmas);
1307   pidObj->SetCompat(compat);
1308   pidObj->SetTPC(kTRUE);
1309   pidObj->SetTOF(kTRUE);
1310   pidObj->SetPCompatTOF(1.5);
1311   pidObj->SetSigmaForTPCCompat(3.);
1312   pidObj->SetSigmaForTOFCompat(3.);
1313   pidObj->SetOldPid(kTRUE);
1314
1315   SetPidHF(pidObj);
1316   SetUsePID(kTRUE);
1317   SetUseDefaultPID(kFALSE);
1318
1319   SetLowPt(kFALSE);
1320
1321   PrintAll();
1322
1323   delete pidObj;
1324   pidObj=NULL;
1325
1326   return;
1327
1328 }
1329
1330
1331 //---------------------------------------------------------------------------
1332 void AliRDHFCutsD0toKpi::SetStandardCutsPP2011_276TeV() {
1333   //
1334   // STANDARD CUTS USED FOR 2011 pp analysis at 2.76TeV
1335   //
1336
1337   SetName("D0toKpiCutsStandard");
1338   SetTitle("Standard Cuts for D0 analysis in pp2011 at 2.76TeV run");
1339
1340   //
1341   // Track cuts
1342   //
1343   AliESDtrackCuts* esdTrackCuts=new AliESDtrackCuts();
1344   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1345   //default
1346   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1347   esdTrackCuts->SetRequireITSRefit(kTRUE);
1348   esdTrackCuts->SetEtaRange(-0.8,0.8);
1349   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1350                                          AliESDtrackCuts::kAny);
1351  // default is kBoth, otherwise kAny
1352   esdTrackCuts->SetMinDCAToVertexXY(0.);
1353   esdTrackCuts->SetPtRange(0.3,1.e10);
1354
1355   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1356   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1357   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1358
1359   AddTrackCuts(esdTrackCuts);
1360
1361
1362   const Int_t nvars=11;
1363   const Int_t nptbins=13;
1364   Float_t ptbins[nptbins+1];
1365   ptbins[0]=0.;
1366   ptbins[1]=0.5;
1367   ptbins[2]=1.;
1368   ptbins[3]=2.;
1369   ptbins[4]=3.;
1370   ptbins[5]=4.;
1371   ptbins[6]=5.;
1372   ptbins[7]=6.;
1373   ptbins[8]=8.;
1374   ptbins[9]=12.;
1375   ptbins[10]=16.;
1376   ptbins[11]=20.;
1377   ptbins[12]=24.;
1378   ptbins[13]=9999.;
1379
1380   SetPtBins(nptbins+1,ptbins);
1381
1382
1383   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*/
1384                                                   {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*/
1385                                                   {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 */
1386                                                   {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 */
1387                                                   {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 */
1388                                                   {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 */
1389                                                   {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 */
1390                                                   {0.400,0.09,0.85,0.7,0.7,1000.*1E-4,1000.*1E-4,0.,0.85,0.,0.},/* 6<pt<8 */
1391                                                   {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 */
1392                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 12<pt<16 */
1393                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 16<pt<20 */
1394                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.},/* 20<pt<24 */
1395                                                   {0.400,0.09,1.0,0.7,0.7,9999.,9999.,0.,0.,0.,0.}};/* pt>24 */
1396
1397   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1398   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1399   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1400   for (Int_t ibin=0;ibin<nptbins;ibin++){
1401     for (Int_t ivar = 0; ivar<nvars; ivar++){
1402       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1403     }
1404   }
1405   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1406   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1407   delete [] cutsMatrixTransposeStand;
1408
1409
1410   //pid settings
1411   AliAODPidHF* pidObj=new AliAODPidHF();
1412   Int_t mode=1;
1413   const Int_t nlims=2;
1414   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1415   Bool_t compat=kTRUE; //effective only for this mode
1416   Bool_t asym=kTRUE;
1417   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1418   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1419   pidObj->SetMatch(mode);
1420   pidObj->SetPLimit(plims,nlims);
1421   pidObj->SetSigma(sigmas);
1422   pidObj->SetCompat(compat);
1423   pidObj->SetTPC(kTRUE);
1424   pidObj->SetTOF(kTRUE);
1425   pidObj->SetOldPid(kTRUE);
1426   SetPidHF(pidObj);
1427
1428   SetUseDefaultPID(kFALSE); //to use the AliAODPidHF
1429
1430   SetUsePID(kTRUE);
1431
1432   //activate pileup rejection (for pp)
1433   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1434
1435   //Do recalculate the vertex
1436   SetRemoveDaughtersFromPrim(kTRUE);
1437
1438   // Cut on the zvtx
1439   SetMaxVtxZ(10.);
1440
1441   // Use the kFirst selection for tracks with candidate D with pt<2
1442   SetSelectCandTrackSPDFirst(kTRUE,2.);
1443
1444   // Use special cuts for D candidates with pt<2
1445   SetUseSpecialCuts(kTRUE);
1446   SetMaximumPtSpecialCuts(2.);
1447
1448   PrintAll();
1449
1450   delete pidObj;
1451   pidObj=NULL;
1452
1453   return;
1454 }
1455
1456
1457 //---------------------------------------------------------------------------
1458 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
1459   //
1460   // Final CUTS USED FOR 2010 PbPb analysis of central events
1461   // REMEMBER TO EVENTUALLY SWITCH ON
1462   //          SetUseAOD049(kTRUE);
1463   //
1464
1465   SetName("D0toKpiCutsStandard");
1466   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
1467
1468
1469   // CENTRALITY SELECTION
1470   SetMinCentrality(0.);
1471   SetMaxCentrality(80.);
1472   SetUseCentrality(AliRDHFCuts::kCentV0M);
1473
1474
1475
1476   // EVENT CUTS
1477   SetMinVtxContr(1);
1478   // MAX Z-VERTEX CUT
1479   SetMaxVtxZ(10.);
1480   // PILE UP
1481   SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1482   // PILE UP REJECTION
1483   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1484
1485   // TRACKS ON SINGLE TRACKS
1486   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1487   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1488   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1489   esdTrackCuts->SetRequireITSRefit(kTRUE);
1490   //  esdTrackCuts->SetMinNClustersITS(4);
1491   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1492   esdTrackCuts->SetMinDCAToVertexXY(0.);
1493   esdTrackCuts->SetEtaRange(-0.8,0.8);
1494   esdTrackCuts->SetPtRange(0.7,1.e10);
1495
1496   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1497   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1498   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0075*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1499
1500
1501   AddTrackCuts(esdTrackCuts);
1502   // SPD k FIRST for Pt<3 GeV/c
1503   SetSelectCandTrackSPDFirst(kTRUE, 3);
1504
1505   // CANDIDATE CUTS
1506   const Int_t nptbins =13;
1507   const Double_t ptmax = 9999.;
1508   const Int_t nvars=11;
1509   Float_t ptbins[nptbins+1];
1510   ptbins[0]=0.;
1511   ptbins[1]=0.5;
1512   ptbins[2]=1.;
1513   ptbins[3]=2.;
1514   ptbins[4]=3.;
1515   ptbins[5]=4.;
1516   ptbins[6]=5.;
1517   ptbins[7]=6.;
1518   ptbins[8]=8.;
1519   ptbins[9]=12.;
1520   ptbins[10]=16.;
1521   ptbins[11]=20.;
1522   ptbins[12]=24.;
1523   ptbins[13]=ptmax;
1524
1525   SetGlobalIndex(nvars,nptbins);
1526   SetPtBins(nptbins+1,ptbins);
1527   SetMinPtCandidate(2.);
1528
1529   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*/
1530                                                   {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*/
1531                                                   {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 */
1532                                                   {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 */
1533                                                   {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 */
1534                                                   {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 */
1535                                                   {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 */
1536                                                   {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 */
1537                                                   {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 */
1538                                                   {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 */
1539                                                   {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 */
1540                                                   {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 */
1541                                                   {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 */
1542
1543
1544   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1545   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1546   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1547
1548   for (Int_t ibin=0;ibin<nptbins;ibin++){
1549     for (Int_t ivar = 0; ivar<nvars; ivar++){
1550       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1551     }
1552   }
1553
1554   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1555   SetUseSpecialCuts(kTRUE);
1556   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1557   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1558   delete [] cutsMatrixTransposeStand;
1559   cutsMatrixTransposeStand=NULL;
1560
1561   // PID SETTINGS
1562   AliAODPidHF* pidObj=new AliAODPidHF();
1563   //pidObj->SetName("pid4D0");
1564   Int_t mode=1;
1565   const Int_t nlims=2;
1566   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1567   Bool_t compat=kTRUE; //effective only for this mode
1568   Bool_t asym=kTRUE;
1569   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1570   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1571   pidObj->SetMatch(mode);
1572   pidObj->SetPLimit(plims,nlims);
1573   pidObj->SetSigma(sigmas);
1574   pidObj->SetCompat(compat);
1575   pidObj->SetTPC(kTRUE);
1576   pidObj->SetTOF(kTRUE);
1577   pidObj->SetPCompatTOF(2.);
1578   pidObj->SetSigmaForTPCCompat(3.);
1579   pidObj->SetSigmaForTOFCompat(3.);
1580   pidObj->SetOldPid(kTRUE);
1581
1582
1583   SetPidHF(pidObj);
1584   SetUsePID(kTRUE);
1585   SetUseDefaultPID(kFALSE);
1586
1587
1588   PrintAll();
1589
1590
1591   delete pidObj;
1592   pidObj=NULL;
1593
1594   return;
1595
1596 }
1597
1598 //---------------------------------------------------------------------------
1599 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010Peripherals() {
1600   // CUTS USED FOR D0 RAA for the CENTRALITY RANGE 20-80%
1601
1602
1603   SetName("D0toKpiCutsStandard");
1604   SetTitle("Standard Cuts for D0 analysis in PbPb2010 run, for peripheral events");
1605
1606
1607   // CENTRALITY SELECTION
1608   SetMinCentrality(40.);
1609   SetMaxCentrality(80.);
1610   SetUseCentrality(AliRDHFCuts::kCentV0M);
1611
1612   // EVENT CUTS
1613   SetMinVtxContr(1);
1614   // MAX Z-VERTEX CUT
1615   SetMaxVtxZ(10.);
1616   // PILE UP
1617   SetOptPileup(AliRDHFCuts::kNoPileupSelection);
1618   // PILE UP REJECTION
1619   //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
1620
1621   // TRACKS ON SINGLE TRACKS
1622   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1623   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1624   esdTrackCuts->SetRequireTPCRefit(kTRUE);
1625   esdTrackCuts->SetRequireITSRefit(kTRUE);
1626   //  esdTrackCuts->SetMinNClustersITS(4);
1627   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
1628   esdTrackCuts->SetMinDCAToVertexXY(0.);
1629   esdTrackCuts->SetEtaRange(-0.8,0.8);
1630   esdTrackCuts->SetPtRange(0.5,1.e10);
1631
1632   esdTrackCuts->SetMaxDCAToVertexXY(1.);
1633   esdTrackCuts->SetMaxDCAToVertexZ(1.);
1634   esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0025*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
1635
1636
1637   AddTrackCuts(esdTrackCuts);
1638   // SPD k FIRST for Pt<3 GeV/c
1639   SetSelectCandTrackSPDFirst(kTRUE, 3);
1640
1641   // CANDIDATE CUTS
1642   const Int_t nptbins =13;
1643   const Double_t ptmax = 9999.;
1644   const Int_t nvars=11;
1645   Float_t ptbins[nptbins+1];
1646   ptbins[0]=0.;
1647   ptbins[1]=0.5;
1648   ptbins[2]=1.;
1649   ptbins[3]=2.;
1650   ptbins[4]=3.;
1651   ptbins[5]=4.;
1652   ptbins[6]=5.;
1653   ptbins[7]=6.;
1654   ptbins[8]=8.;
1655   ptbins[9]=12.;
1656   ptbins[10]=16.;
1657   ptbins[11]=20.;
1658   ptbins[12]=24.;
1659   ptbins[13]=ptmax;
1660
1661   SetGlobalIndex(nvars,nptbins);
1662   SetPtBins(nptbins+1,ptbins);
1663   SetMinPtCandidate(0.);
1664
1665   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*/
1666                                                   {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*/
1667                                                   {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 */
1668                                                   {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 */
1669                                                   {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 */
1670                                                   {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 */
1671                                                   {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 */
1672                                                   {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 */
1673                                                   {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 */
1674                                                   {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 */
1675                                                   {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 */
1676                                                   {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 */
1677                                                   {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 */
1678
1679
1680   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
1681   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
1682   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
1683
1684   for (Int_t ibin=0;ibin<nptbins;ibin++){
1685     for (Int_t ivar = 0; ivar<nvars; ivar++){
1686       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
1687     }
1688   }
1689
1690   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
1691   SetUseSpecialCuts(kTRUE);
1692   SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
1693   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
1694   delete [] cutsMatrixTransposeStand;
1695   cutsMatrixTransposeStand=NULL;
1696
1697   // PID SETTINGS
1698   AliAODPidHF* pidObj=new AliAODPidHF();
1699   //pidObj->SetName("pid4D0");
1700   Int_t mode=1;
1701   const Int_t nlims=2;
1702   Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
1703   Bool_t compat=kTRUE; //effective only for this mode
1704   Bool_t asym=kTRUE;
1705   Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
1706   pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
1707   pidObj->SetMatch(mode);
1708   pidObj->SetPLimit(plims,nlims);
1709   pidObj->SetSigma(sigmas);
1710   pidObj->SetCompat(compat);
1711   pidObj->SetTPC(kTRUE);
1712   pidObj->SetTOF(kTRUE);
1713   pidObj->SetPCompatTOF(2.);
1714   pidObj->SetSigmaForTPCCompat(3.);
1715   pidObj->SetSigmaForTOFCompat(3.);
1716   pidObj->SetOldPid(kTRUE);
1717
1718   SetPidHF(pidObj);
1719   SetUsePID(kTRUE);
1720   SetUseDefaultPID(kFALSE);
1721
1722   SetLowPt(kTRUE,2.);
1723
1724   PrintAll();
1725
1726
1727   delete pidObj;
1728   pidObj=NULL;
1729
1730   return;
1731
1732 }
1733
1734
1735 //---------------------------------------------------------------------------
1736 void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2011() {
1737
1738   // Default 2010 PbPb cut object
1739   SetStandardCutsPbPb2010();
1740   AliAODPidHF *pidobj=GetPidHF();
1741
1742   pidobj->SetOldPid(kFALSE);
1743
1744   //
1745   // Enable all 2011 PbPb run triggers
1746   //
1747   SetTriggerClass("");
1748   ResetMaskAndEnableMBTrigger();
1749   EnableCentralTrigger();
1750   EnableSemiCentralTrigger();
1751
1752 }
1753
1754 //---------------------------------------------------------------------------
1755 Int_t AliRDHFCutsD0toKpi::IsSelectedCombPID(AliAODRecoDecayHF* d)
1756 {
1757   // ############################################################
1758   //
1759   // Apply Bayesian PID selection
1760   // Implementation of Bayesian PID: Jeremy Wilkinson
1761   //
1762   // ############################################################
1763
1764   if (!fUsePID || !d) return 3;
1765
1766
1767   if (fBayesianStrategy == kBayesWeightNoFilter) {
1768      //WeightNoFilter: Accept all particles (no PID cut) but fill mass histos with weights in task
1769      CalculateBayesianWeights(d);
1770      return 3;
1771
1772 }
1773
1774
1775
1776   Int_t isD0 = 0;
1777   Int_t isD0bar = 0;
1778   Int_t returnvalue = 0;
1779
1780   Int_t isPosKaon = 0, isNegKaon = 0;
1781
1782   //Bayesian methods used here check for ID of kaon, and whether it is positive or negative.
1783
1784
1785   Bool_t checkPIDInfo[2] = {kTRUE, kTRUE};
1786   AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1787
1788   if ((aodtrack[0]->Charge()*aodtrack[1]->Charge()) != -1) {
1789     return 0;  //reject if charges not opposite
1790   }
1791    Double_t momentumpositive=0., momentumnegative=0.;   //Used in "prob > prior" method
1792   for (Int_t daught = 0; daught < 2; daught++) {
1793     //Loop over prongs
1794
1795       if (aodtrack[daught]->Charge() == -1) {
1796          momentumnegative = aodtrack[daught]->P();
1797       }
1798       if (aodtrack[daught]->Charge() == +1) {
1799          momentumpositive = aodtrack[daught]->P();
1800       }
1801           if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1802              checkPIDInfo[daught] = kFALSE;
1803              continue;
1804           }
1805
1806    }
1807
1808    //Loop over daughters ends here
1809
1810    if (!checkPIDInfo[0] && !checkPIDInfo[1]) {
1811       return 0;      //Reject if both daughters lack both TPC and TOF info
1812    }
1813
1814
1815    CalculateBayesianWeights(d);        //Calculates all Bayesian probabilities for both positive and negative tracks
1816     //      Double_t prob0[AliPID::kSPECIES];
1817     //      fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), prob0);
1818    ///Three possible Bayesian probability cuts: Picked using SetBayesianCondition(int).
1819    switch (fBayesianCondition) {
1820       ///A: Standard max. probability method (accept most likely species)
1821       case kMaxProb:
1822                         if (TMath::MaxElement(AliPID::kSPECIES, fWeightsPositive) == fWeightsPositive[AliPID::kKaon]) { //If highest probability lies with kaon
1823                           isPosKaon = 1;  //flag [daught] as a kaon
1824                         }
1825
1826                         if (TMath::MaxElement(AliPID::kSPECIES, fWeightsNegative) == fWeightsNegative[AliPID::kKaon]) { //If highest probability lies with kaon
1827                           isNegKaon = 1;  //flag [daught] as a kaon
1828                         }
1829                     break;
1830       ///B: Accept if probability greater than prior
1831       case kAbovePrior:
1832
1833                         if (fWeightsNegative[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1834                                                                         GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumnegative)))) {  //Retrieves relevant prior, gets value at momentum
1835                         isNegKaon = 1;
1836                         }
1837                         if (fWeightsPositive[AliPID::kKaon] > (fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->
1838                                                                         GetBinContent(fPidHF->GetPidCombined()->GetPriorDistribution(AliPID::kKaon)->FindBin(momentumpositive)))) {  //Retrieves relevant prior, gets value at momentum
1839                         isPosKaon = 1;
1840                         }
1841
1842                 break;
1843
1844       ///C: Accept if probability greater than user-defined threshold
1845       case kThreshold:
1846          if (fWeightsNegative[AliPID::kKaon] > fProbThreshold) {
1847             isNegKaon = 1;
1848          }
1849          if (fWeightsPositive[AliPID::kKaon] > fProbThreshold) {
1850             isPosKaon = 1;
1851          }
1852
1853          break;
1854   }
1855
1856
1857      //Momentum-based selection (also applied to filtered weighted method)
1858
1859      if (fBayesianStrategy == kBayesMomentum || fBayesianCondition == kBayesWeight) {
1860          if (isNegKaon && isPosKaon) { // If both are kaons, reject
1861             isD0 = 0;
1862             isD0bar = 0;
1863          } else if (isNegKaon) {       //If negative kaon present, D0
1864             isD0 = 1;
1865          } else if (isPosKaon) {       //If positive kaon present, D0bar
1866             isD0bar = 1;
1867          } else {                      //If neither ID'd as kaon, subject to extra tests
1868             isD0 = 1;
1869             isD0bar = 1;
1870             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
1871                if (aodtrack[0]->Charge() == -1) {
1872             isD0 = 0;  //Check charge and reject based on pion hypothesis
1873                }
1874                if (aodtrack[0]->Charge() == 1) {
1875             isD0bar = 0;
1876                }
1877             }
1878             if (aodtrack[1]->P() < 1.5 && (fPidHF->CheckStatus(aodtrack[1], "TOF"))) {
1879                if (aodtrack[1]->Charge() == -1) {
1880             isD0 = 0;
1881                }
1882                if (aodtrack[1]->Charge() == 1) {
1883             isD0bar = 0;
1884                }
1885             }
1886          }
1887
1888
1889
1890          if (isD0 && isD0bar) {
1891             returnvalue = 3;
1892          }
1893          if (isD0 && !isD0bar) {
1894             returnvalue = 1;
1895          }
1896          if (!isD0 && isD0bar) {
1897             returnvalue = 2;
1898          }
1899          if (!isD0 && !isD0bar) {
1900             returnvalue = 0;
1901          }
1902      }
1903
1904     //Simple Bayesian
1905     if (fBayesianStrategy == kBayesSimple) {
1906
1907          if (isPosKaon && isNegKaon)   {  //If both are ID'd as kaons, accept as possible
1908                returnvalue = 3;
1909             } else if (isNegKaon)   {     //If negative kaon, D0
1910                returnvalue = 1;
1911             } else if (isPosKaon)   {     //If positive kaon, D0-bar
1912                returnvalue = 2;
1913             } else {
1914                returnvalue = 0;  //If neither kaon, reject
1915             }
1916     }
1917
1918   return returnvalue;
1919
1920
1921
1922 }
1923
1924
1925
1926 //---------------------------------------------------------------------------
1927 void AliRDHFCutsD0toKpi::CalculateBayesianWeights(AliAODRecoDecayHF* d)
1928
1929 {
1930   //Function to compute weights for Bayesian method
1931
1932   AliAODTrack *aodtrack[2] = {(AliAODTrack*)d->GetDaughter(0), (AliAODTrack*)d->GetDaughter(1)};
1933   if ((aodtrack[0]->Charge() * aodtrack[1]->Charge()) != -1) {
1934     return;  //Reject if charges do not oppose
1935   }
1936   for (Int_t daught = 0; daught < 2; daught++) {
1937     //Loop over prongs
1938
1939     if (!(fPidHF->CheckStatus(aodtrack[daught], "TPC")) && !(fPidHF->CheckStatus(aodtrack[daught], "TOF"))) {
1940       continue;
1941     }
1942
1943
1944     // identify kaon, define weights
1945     if (aodtrack[daught]->Charge() == +1) {
1946     fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsPositive);
1947     }
1948
1949     if (aodtrack[daught]->Charge() == -1) {
1950     fPidHF->GetPidCombined()->ComputeProbabilities(aodtrack[daught], fPidHF->GetPidResponse(), fWeightsNegative);
1951     }
1952   }
1953 }
1954
1955
1956