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