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