PWGHF/vertexingHF/AliRDHFCutsD0toKpi.cxx

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