PWG3/vertexingHF/AliRDHFCutsDStartoKpipi.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 /////////////////////////////////////////////////////////////
  19 //
  20 // Class for cuts on AOD reconstructed DStar->Kpipi
  21 //
  22 // Author: A.Grelli, alessandro.grelli@uu.nl
  23 //
  24 // PID method implemented by   Y.Wang, yifei@physi.uni-heidelberg.de
  25 //
  26 /////////////////////////////////////////////////////////////
  27
  28 #include <TDatabasePDG.h>
  29 #include <Riostream.h>
  30
  31 #include "AliAODRecoDecayHF2Prong.h"
  32 #include "AliAODRecoCascadeHF.h"
  33 #include "AliRDHFCutsD0toKpi.h"
  34 #include "AliRDHFCutsDStartoKpipi.h"
  35 #include "AliAODTrack.h"
  36 #include "AliESDtrack.h"
  37 #include "AliAODPid.h"
  38 #include "AliTPCPIDResponse.h"
  39 #include "AliAODVertex.h"
  40 #include "AliESDVertex.h"
  41
  42 ClassImp(AliRDHFCutsDStartoKpipi)
  43
  44 //--------------------------------------------------------------------------
  45 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const char* name) :
  46   AliRDHFCuts(name),
  47   fTrackCutsSoftPi(0)
  48 {
  49   //
  50   // Default Constructor
  51   //
  52
  53   Int_t nvars=16;
  54   SetNVars(nvars);
  55   TString varNames[16]={
  56     "inv. mass [GeV]",
  57     "dca [cm]",
  58     "cosThetaStar",
  59     "pTK [GeV/c]",
  60     "pTPi [GeV/c]",
  61     "d0K [cm]",
  62     "d0Pi [cm]",
  63     "d0d0 [cm^2]",
  64     "cosThetaPoint",
  65     "inv. mass half width of D* [GeV]",
  66     "half width of (M_Kpipi-M_D0) [GeV]",
  67     "PtMin of pi_s [GeV/c]",
  68     "PtMax of pi_s [GeV/c]",
  69     "theta, angle between the pi_s and decay plane of the D0 [rad]",
  70     "|cosThetaPointXY|",
  71     "NormDecayLenghtXY"};
  72   Bool_t isUpperCut[16]={
  73     kTRUE,
  74     kTRUE,
  75     kTRUE,
  76     kFALSE,
  77     kFALSE,
  78     kTRUE,
  79     kTRUE,
  80     kTRUE,
  81     kFALSE,
  82     kTRUE,
  83     kTRUE,
  84     kTRUE,
  85     kTRUE,
  86     kFALSE,
  87     kFALSE,
  88     kFALSE};
  89   SetVarNames(nvars,varNames,isUpperCut);
  90   Bool_t forOpt[16]={
  91     kFALSE,
  92     kTRUE,
  93     kTRUE,
  94     kFALSE,
  95     kFALSE,
  96     kFALSE,
  97     kFALSE,
  98     kTRUE,
  99     kTRUE,
 100     kFALSE,
 101     kTRUE,
 102     kFALSE,
 103     kFALSE,
 104     kFALSE,
 105     kFALSE,
 106     kFALSE};
 107   SetVarsForOpt(5,forOpt);
 108   Float_t limits[2]={0,999999999.};
 109   SetPtBins(2,limits);
 110 }
 111 //--------------------------------------------------------------------------
 112 AliRDHFCutsDStartoKpipi::AliRDHFCutsDStartoKpipi(const AliRDHFCutsDStartoKpipi &source) :
 113   AliRDHFCuts(source),
 114   fTrackCutsSoftPi(0)
 115 {
 116   //
 117   // Copy constructor
 118   //
 119
 120   if(source.GetTrackCutsSoftPi()) AddTrackCutsSoftPi(source.GetTrackCutsSoftPi());
 121
 122 }
 123 //--------------------------------------------------------------------------
 124 AliRDHFCutsDStartoKpipi &AliRDHFCutsDStartoKpipi::operator=(const AliRDHFCutsDStartoKpipi &source)
 125 {
 126   //
 127   // assignment operator
 128   //
 129   if(&source == this) return *this;
 130
 131   AliRDHFCuts::operator=(source);
 132   if(source.GetTrackCutsSoftPi()) AddTrackCutsSoftPi(source.GetTrackCutsSoftPi());
 133
 134   return *this;
 135 }
 136
 137
 138 //---------------------------------------------------------------------------
 139 void AliRDHFCutsDStartoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
 140   //
 141   // Fills in vars the values of the variables
 142   //
 143   if(nvars!=fnVarsForOpt) {
 144     printf("AliRDHFCutsDStartoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
 145     return;
 146   }
 147
 148
 149   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)d;
 150
 151   AliAODTrack *softPi = (AliAODTrack*)dstarD0pi->GetBachelor();
 152
 153   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 154
 155    Int_t iter=-1;
 156   if(fVarsForOpt[0]){
 157     iter++;
 158     if(TMath::Abs(pdgdaughters[0])==211) {
 159       vars[iter]=dd->InvMassD0();
 160     } else {
 161       vars[iter]=dd->InvMassD0bar();
 162     }
 163   }
 164   if(fVarsForOpt[1]){
 165     iter++;
 166     vars[iter]=dd->GetDCA();
 167   }
 168   if(fVarsForOpt[2]){
 169     iter++;
 170     if(TMath::Abs(pdgdaughters[0])==211) {
 171       vars[iter] = dd->CosThetaStarD0();
 172     } else {
 173       vars[iter] = dd->CosThetaStarD0bar();
 174     }
 175   }
 176   if(fVarsForOpt[3]){
 177     iter++;
 178    if(TMath::Abs(pdgdaughters[0])==321) {
 179      vars[iter]=dd->PtProng(0);
 180    }
 181    else{
 182      vars[iter]=dd->PtProng(1);
 183    }
 184   }
 185   if(fVarsForOpt[4]){
 186     iter++;
 187    if(TMath::Abs(pdgdaughters[0])==211) {
 188      vars[iter]=dd->PtProng(0);
 189    }
 190    else{
 191      vars[iter]=dd->PtProng(1);
 192    }
 193   }
 194   if(fVarsForOpt[5]){
 195     iter++;
 196     if(TMath::Abs(pdgdaughters[0])==321) {
 197      vars[iter]=dd->Getd0Prong(0);
 198    }
 199    else{
 200      vars[iter]=dd->Getd0Prong(1);
 201    }
 202   }
 203   if(fVarsForOpt[6]){
 204     iter++;
 205      if(TMath::Abs(pdgdaughters[0])==211) {
 206      vars[iter]=dd->Getd0Prong(0);
 207    }
 208    else{
 209      vars[iter]=dd->Getd0Prong(1);
 210    }
 211   }
 212   if(fVarsForOpt[7]){
 213     iter++;
 214     vars[iter]= dd->Prodd0d0();
 215   }
 216   if(fVarsForOpt[8]){
 217     iter++;
 218     vars[iter]=dd->CosPointingAngle();
 219   }
 220   if(fVarsForOpt[9]){
 221     iter++;
 222     vars[iter]=dstarD0pi->InvMassDstarKpipi();
 223   }
 224   if(fVarsForOpt[10]){
 225     iter++;
 226     vars[iter]=dstarD0pi->DeltaInvMass();
 227   }
 228   if(fVarsForOpt[11]){
 229     iter++;
 230     vars[iter] = softPi->Pt();
 231   }
 232   if(fVarsForOpt[12]){
 233     iter++;
 234     vars[iter] = softPi->Pt();
 235   }
 236   if(fVarsForOpt[13]){
 237     iter++;
 238     vars[iter] =dstarD0pi->AngleD0dkpPisoft();
 239   }
 240   if(fVarsForOpt[14]){
 241     iter++;
 242     vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
 243   }
 244   if(fVarsForOpt[15]){
 245     iter++;
 246     vars[iter]=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
 247   }
 248
 249   return;
 250 }
 251 //---------------------------------------------------------------------------
 252 Int_t AliRDHFCutsDStartoKpipi::IsSelected(TObject* obj,Int_t selectionLevel) {
 253   //
 254   // Apply selection for D*.
 255   //
 256   if(!fCutsRD){
 257     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 258     return 0;
 259   }
 260
 261   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
 262   if(!d){
 263     cout<<"AliAODRecoCascadeHF null"<<endl;
 264     return 0;
 265   }
 266
 267   Double_t ptD=d->Pt();
 268   if(ptD<fMinPtCand) return 0;
 269   if(ptD>fMaxPtCand) return 0;
 270
 271
 272   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)d->Get2Prong();
 273   if(!dd){
 274     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 275     return 0;
 276   }
 277
 278   if(dd->HasBadDaughters()) return 0;
 279
 280   AliAODTrack *b = (AliAODTrack*)d->GetBachelor();
 281   if(fTrackCutsSoftPi && fTrackCutsSoftPi->GetRequireTPCRefit()){
 282     if(!(b->TestFilterMask(BIT(4)))) return 0;
 283   }
 284
 285   Int_t returnvalue=1;
 286   Int_t returnvaluePID=3;
 287
 288
 289   // selection on candidate
 290   if(selectionLevel==AliRDHFCuts::kAll ||
 291      selectionLevel==AliRDHFCuts::kCandidate) {
 292
 293     Double_t pt=d->Pt();
 294     Int_t ptbin=PtBin(pt);
 295
 296     // DStarMass and D0mass
 297     Double_t mDSPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 298     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 299     // delta mass PDG
 300     Double_t deltaPDG = mDSPDG-mD0PDG;
 301
 302     // Half width DStar mass
 303     if(TMath::Abs(mDSPDG - (d->InvMassDstarKpipi()))>fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
 304     // Half width Delta mass
 305
 306     if(TMath::Abs(deltaPDG-(d->DeltaInvMass())) > fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
 307
 308     // cut on soft pion pt
 309     if(b->Pt() < fCutsRD[GetGlobalIndex(11,ptbin)] || b->Pt() > fCutsRD[GetGlobalIndex(12,ptbin)]) return 0;
 310     // cut on the angle between D0 decay plane and soft pion
 311     if(d->AngleD0dkpPisoft() > fCutsRD[GetGlobalIndex(13,ptbin)]) return 0;
 312
 313     // select D0 that passes D* cuts
 314     returnvalue = IsD0FromDStarSelected(pt,dd,selectionLevel);
 315     if((b->Charge()==+1 && returnvalue==2) || (b->Charge()==-1 && returnvalue==1)) return 0;
 316
 317   }
 318
 319   // selection on PID
 320   if(selectionLevel==AliRDHFCuts::kAll ||
 321      selectionLevel==AliRDHFCuts::kCandidate ||
 322      selectionLevel==AliRDHFCuts::kPID) {
 323     returnvaluePID = IsSelectedPID(d);
 324   }
 325   if(returnvaluePID!=3) returnvalue =0;
 326
 327
 328   // selection on daughter tracks
 329   if(selectionLevel==AliRDHFCuts::kAll ||
 330      selectionLevel==AliRDHFCuts::kTracks) {
 331     if(!AreDaughtersSelected(dd)) return 0;
 332     if(fTrackCutsSoftPi) {
 333       AliAODVertex *vAOD = d->GetPrimaryVtx();
 334       Double_t pos[3],cov[6];
 335       vAOD->GetXYZ(pos);
 336       vAOD->GetCovarianceMatrix(cov);
 337       const AliESDVertex vESD(pos,cov,100.,100);
 338       if(!IsDaughterSelected(b,&vESD,fTrackCutsSoftPi)) return 0;
 339     }
 340   }
 341
 342   return returnvalue;
 343
 344 }
 345 //_________________________________________________________________________________________________
 346 Int_t AliRDHFCutsDStartoKpipi::IsD0FromDStarSelected(Double_t pt, TObject* obj,Int_t selectionLevel) const {
 347   //
 348   // Apply selection for D0 from D*. The selection in on D0 prongs
 349   //
 350
 351   if(!fCutsRD){
 352     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 353     return 0;
 354   }
 355
 356   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)obj;
 357
 358   if(!dd){
 359     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 360     return 0;
 361   }
 362
 363   // selection on daughter tracks is done in IsSelected()
 364
 365   Int_t returnvalue=1;
 366
 367   // selection on candidate
 368   if(selectionLevel==AliRDHFCuts::kAll ||
 369      selectionLevel==AliRDHFCuts::kCandidate) {
 370
 371     // D0 mass
 372     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 373     // delta mass PDG
 374
 375     Int_t ptbin=PtBin(pt);
 376
 377     Double_t mD0,mD0bar,ctsD0,ctsD0bar;
 378
 379     Int_t okD0     =0;
 380     Int_t okD0bar  =0;
 381     okD0=1; okD0bar=1;
 382
 383     if(dd->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
 384     if(dd->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
 385
 386     if(!okD0 && !okD0bar) return 0;
 387
 388     if(TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
 389        TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
 390     if(TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
 391        TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
 392     if(!okD0 && !okD0bar) return 0;
 393
 394     if(dd->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
 395
 396     dd->InvMassD0(mD0,mD0bar);
 397     if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
 398     if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
 399     if(!okD0 && !okD0bar) return 0;
 400
 401     dd->CosThetaStarD0(ctsD0,ctsD0bar);
 402     if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
 403     if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
 404     if(!okD0 && !okD0bar) return 0;
 405
 406     if(dd->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
 407
 408     if(dd->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
 409
 410     if(TMath::Abs(dd->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(14,ptbin)])  return 0;
 411
 412     Double_t normalDecayLengXY=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
 413     if (normalDecayLengXY < fCutsRD[GetGlobalIndex(15, ptbin)]) return 0;
 414
 415     if (okD0) returnvalue=1; //cuts passed as D0
 416     if (okD0bar) returnvalue=2; //cuts passed as D0bar
 417     if (okD0 && okD0bar) returnvalue=3; //both
 418   }
 419
 420   return returnvalue;
 421 }
 422 //----------------------------------------------------------------------------------
 423 Bool_t AliRDHFCutsDStartoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 424 {
 425   //
 426   // D* fiducial acceptance region
 427   //
 428
 429   if(pt > 5.) {
 430     // applying cut for pt > 5 GeV
 431     AliDebug(4,Form("pt of D* = %f (> 5), cutting at |y| < 0.8\n",pt));
 432     if (TMath::Abs(y) > 0.8){
 433       return kFALSE;
 434     }
 435   } else {
 436     // appliying smooth cut for pt < 5 GeV
 437     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 438     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 439     AliDebug(2,Form("pt of D* = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 440     if (y < minFiducialY || y > maxFiducialY){
 441       return kFALSE;
 442     }
 443   }
 444
 445   return kTRUE;
 446 }
 447
 448 //_______________________________________________________________________________-
 449 Int_t AliRDHFCutsDStartoKpipi::IsSelectedPID(AliAODRecoDecayHF* obj)
 450 {
 451   //
 452   // PID method, n signa approach default
 453   //
 454
 455   if(!fUsePID) return 3;
 456
 457   AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)obj;
 458   if(!dstar){
 459     cout<<"AliAODRecoCascadeHF null"<<endl;
 460     return 0;
 461   }
 462   AliAODRecoDecayHF2Prong* d0 = (AliAODRecoDecayHF2Prong*)dstar->Get2Prong();
 463   if(!d0){
 464     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 465     return 0;
 466   }
 467
 468   //  here the PID
 469   AliAODTrack *pos = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(0);
 470   AliAODTrack *neg = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(1);
 471
 472   if (dstar->Charge()>0){
 473     if(!SelectPID(pos,2)) return 0;//pion+
 474     if(!SelectPID(neg,3)) return 0;//kaon-
 475   }else{
 476     if(!SelectPID(pos,3)) return 0;//kaon+
 477     if(!SelectPID(neg,2)) return 0;//pion-
 478   }
 479
 480   return 3;
 481 }
 482
 483 //_______________________________________________________________________________-
 484 Int_t AliRDHFCutsDStartoKpipi::SelectPID(AliAODTrack *track, Int_t type)
 485 {
 486   //
 487   //  here the PID
 488
 489   Bool_t isParticle=kTRUE;
 490
 491   if(fPidHF->GetMatch()==1){//n-sigma
 492     Bool_t TPCon=TMath::Abs(2)>1e-4?kTRUE:kFALSE;
 493     Bool_t TOFon=TMath::Abs(3)>1e-4?kTRUE:kFALSE;
 494
 495     Bool_t isTPC=kTRUE;
 496     Bool_t isTOF=kTRUE;
 497
 498     if (TPCon){//TPC
 499       if(fPidHF->CheckStatus(track,"TPC")){
 500         if(type==2) isTPC=fPidHF->IsPionRaw(track,"TPC");
 501         if(type==3) isTPC=fPidHF->IsKaonRaw(track,"TPC");
 502       }
 503     }
 504     if (TOFon){//TOF
 505       if(fPidHF->CheckStatus(track,"TOF")){
 506         if(type==2) isTOF=fPidHF->IsPionRaw(track,"TOF");
 507         if(type==3) isTOF=fPidHF->IsKaonRaw(track,"TOF");
 508       }
 509     }
 510     isParticle = isTPC&&isTOF;
 511   }
 512
 513   if(fPidHF->GetMatch()==2){//bayesian
 514     //Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 515     Double_t prob[5]={1.,1.,1.,1.,1.};
 516
 517     //fPidHF->SetPriors(priors);
 518     fPidHF->BayesianProbability(track,prob);
 519
 520     Double_t max=0.;
 521     Int_t k=-1;
 522     for (Int_t i=0; i<5; i++) {
 523       if (prob[i]>max) {k=i; max=prob[i];}
 524     }
 525     isParticle = Bool_t(k==type);
 526   }
 527
 528   return isParticle;
 529
 530 }
 531 //__________________________________________________________________________________-
 532 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPP2010() {
 533   //
 534   //STANDARD CUTS USED FOR 2010 pp analysis
 535   //
 536   // Need to be updated for the final cut version
 537   //
 538
 539   SetName("DStartoD0piCutsStandard");
 540   SetTitle("Standard Cuts for D* analysis");
 541
 542   // PILE UP REJECTION
 543   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 544
 545   // EVENT CUTS
 546   SetMinVtxContr(1);
 547
 548   // CUTS ON SINGLE TRACKS
 549   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 550   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 551   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 552   esdTrackCuts->SetRequireITSRefit(kTRUE);
 553   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 554   esdTrackCuts->SetMinDCAToVertexXY(0.);
 555   esdTrackCuts->SetEtaRange(-0.8,0.8);
 556   esdTrackCuts->SetPtRange(0.3,1.e10);
 557
 558   // CUTS on SOFT PION
 559   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
 560   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
 561   esdSoftPicuts->SetRequireTPCRefit(kFALSE);
 562   esdSoftPicuts->SetRequireITSRefit(kFALSE);
 563   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 564                                           AliESDtrackCuts::kAny);
 565   esdSoftPicuts->SetPtRange(0.0,1.e10);
 566
 567   AddTrackCuts(esdTrackCuts);
 568   AddTrackCutsSoftPi(esdSoftPicuts);
 569
 570   const Int_t nptbins =13;
 571   const Double_t ptmax = 9999.;
 572   const Int_t nvars=16;
 573   Float_t ptbins[nptbins+1];
 574   ptbins[0]=0.;
 575   ptbins[1]=0.5;
 576   ptbins[2]=1.;
 577   ptbins[3]=2.;
 578   ptbins[4]=3.;
 579   ptbins[5]=4.;
 580   ptbins[6]=5.;
 581   ptbins[7]=6.;
 582   ptbins[8]=7.;
 583   ptbins[9]=8.;
 584   ptbins[10]=12.;
 585   ptbins[11]=16.;
 586   ptbins[12]=24.;
 587   ptbins[13]=ptmax;
 588
 589   SetGlobalIndex(nvars,nptbins);
 590   SetPtBins(nptbins+1,ptbins);
 591
 592   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-2000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* pt<0.5*/
 593                                                   {0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-16000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* 0.5<pt<1*/
 594                                                   {0.7,400.*1E-4,0.8,0.7,0.7,400.*1E-4,400.*1E-4,-36000.*1E-8,0.82,0.3,0.1,0.05,100,0.5,-1.,0.},/* 1<pt<2 */
 595                                                   {0.7,200.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-16000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 2<pt<3 */
 596                                                   {0.7,500.*1E-4,0.8,1.0,1.0,420.*1E-4,560.*1E-4,-6500.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 3<pt<4 */
 597                                                   {0.7,800.*1E-4,0.9,1.2,1.2,700.*1E-4,700.*1E-4,1000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 4<pt<5 */
 598                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,800.*1E-4,800.*1E-4,50000.*1E-8,0.8,0.3,0.1,0.05,100,0.5,-1.,0.},/* 5<pt<6 */
 599                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 6<pt<7 */
 600                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 7<pt<8 */
 601                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,600000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 8<pt<12 */
 602                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 12<pt<16 */
 603                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 16<pt<20 */
 604                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.}};/* pt>24 */
 605
 606
 607   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 608   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 609   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 610
 611   for (Int_t ibin=0;ibin<nptbins;ibin++){
 612     for (Int_t ivar = 0; ivar<nvars; ivar++){
 613       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
 614     }
 615   }
 616
 617   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 618
 619   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
 620   delete [] cutsMatrixTransposeStand;
 621   cutsMatrixTransposeStand=NULL;
 622
 623   // PID SETTINGS FOR D* analysis
 624   AliAODPidHF* pidObj=new AliAODPidHF();
 625   //pidObj->SetName("pid4DSatr");
 626   Int_t mode=1;
 627   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 628   pidObj->SetPriors(priors);
 629   pidObj->SetMatch(mode);
 630   pidObj->SetSigma(0,2); // TPC
 631   pidObj->SetSigma(3,3); // TOF
 632   pidObj->SetTPC(kTRUE);
 633   pidObj->SetTOF(kTRUE);
 634
 635   SetPidHF(pidObj);
 636   SetUsePID(kTRUE);
 637
 638   PrintAll();
 639
 640   delete pidObj;
 641   pidObj=NULL;
 642
 643   return;
 644 }
 645 //_____________________________________________________________________________-
 646 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPbPb2010(){
 647   //
 648   // TEMPORARY, WORK IN PROGRESS ... BUT WORKING!
 649   //
 650   //  Lead Lead
 651   //
 652
 653   SetName("DStartoD0piCutsStandard");
 654   SetTitle("Standard Cuts for D* analysis in PbPb 2010");
 655
 656   // EVENT CUTS
 657   SetMinVtxContr(1);
 658
 659   // CUTS ON SINGLE TRACKS
 660   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 661   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 662   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 663   esdTrackCuts->SetRequireITSRefit(kTRUE);
 664   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 665   esdTrackCuts->SetMinDCAToVertexXY(0.);
 666   esdTrackCuts->SetEtaRange(-0.8,0.8);
 667   esdTrackCuts->SetPtRange(0.3,1.e10);
 668
 669   // CUTS on SOFT PION
 670   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
 671   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
 672   esdSoftPicuts->SetRequireTPCRefit(kTRUE);
 673   esdSoftPicuts->SetRequireITSRefit(kTRUE);
 674   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 675                                           AliESDtrackCuts::kAny); //test d0 asimmetry
 676   esdSoftPicuts->SetPtRange(0.25,5);
 677
 678   AddTrackCuts(esdTrackCuts);
 679   AddTrackCutsSoftPi(esdSoftPicuts);
 680
 681   const Int_t nptbins =13;
 682   const Double_t ptmax = 9999.;
 683   const Int_t nvars=16;
 684   Float_t ptbins[nptbins+1];
 685   ptbins[0]=0.;
 686   ptbins[1]=0.5;
 687   ptbins[2]=1.;
 688   ptbins[3]=2.;
 689   ptbins[4]=3.;
 690   ptbins[5]=4.;
 691   ptbins[6]=5.;
 692   ptbins[7]=6.;
 693   ptbins[8]=7.;
 694   ptbins[9]=8.;
 695   ptbins[10]=12.;
 696   ptbins[11]=16.;
 697   ptbins[12]=24.;
 698   ptbins[13]=ptmax;
 699
 700   SetGlobalIndex(nvars,nptbins);
 701   SetPtBins(nptbins+1,ptbins);
 702
 703   Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-2000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* pt<0.5*/
 704                                                   {0.7,220.*1E-4,0.7,0.21,0.21,500.*1E-4,500.*1E-4,-16000.*1E-8,0.85,0.3,0.1,0.05,100,0.5,-1.,0.},/* 0.5<pt<1*/
 705                                                   {0.7,400.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-36000.*1E-8,0.82,0.3,0.1,0.05,100,0.5,-1.,0.},/* 1<pt<2 */
 706                                                   {0.7,200.*1E-4,0.8,0.7,0.7,800.*1E-4,800.*1E-4,-16000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 2<pt<3 */
 707                                                   {0.7,500.*1E-4,0.8,1.0,1.0,420.*1E-4,560.*1E-4,-6500.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 3<pt<4 */
 708                                                   {0.7,800.*1E-4,0.9,1.2,1.2,700.*1E-4,700.*1E-4,1000.*1E-8,0.9,0.3,0.1,0.05,100,0.5,-1.,0.},/* 4<pt<5 */
 709                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,800.*1E-4,800.*1E-4,50000.*1E-8,0.8,0.3,0.1,0.05,100,0.5,-1.,0.},/* 5<pt<6 */
 710                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 6<pt<7 */
 711                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,100000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 7<pt<8 */
 712                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1000.*1E-4,1000.*1E-4,600000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 8<pt<12 */
 713                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 12<pt<16 */
 714                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.},/* 16<pt<24 */
 715                                                   {0.7,1000.*1E-4,1.0,1.0,1.0,1500.*1E-4,1500.*1E-4,1000000.*1E-8,0.7,0.3,0.1,0.05,100,0.5,-1.,0.}};/* pt>24 */
 716
 717
 718   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 719   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 720   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 721
 722   for (Int_t ibin=0;ibin<nptbins;ibin++){
 723     for (Int_t ivar = 0; ivar<nvars; ivar++){
 724       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
 725     }
 726   }
 727
 728   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 729
 730   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
 731   delete [] cutsMatrixTransposeStand;
 732   cutsMatrixTransposeStand=NULL;
 733
 734   // PID SETTINGS
 735   AliAODPidHF* pidObj=new AliAODPidHF();
 736   // pidObj->SetName("pid4DSatr");
 737   Int_t mode=1;
 738   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 739   pidObj->SetPriors(priors);
 740   pidObj->SetMatch(mode);
 741   pidObj->SetSigma(0,2); // TPC
 742   pidObj->SetSigma(3,3); // TOF
 743   pidObj->SetTPC(kTRUE);
 744   pidObj->SetTOF(kTRUE);
 745
 746   SetPidHF(pidObj);
 747   SetUsePID(kTRUE);
 748
 749   PrintAll();
 750
 751   delete pidObj;
 752   pidObj=NULL;
 753
 754   return;
 755
 756 }