PWGHF/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()) {
 133     delete fTrackCutsSoftPi;
 134     fTrackCutsSoftPi = new AliESDtrackCuts(*(source.GetTrackCutsSoftPi()));
 135   }
 136
 137   return *this;
 138 }
 139
 140
 141 //---------------------------------------------------------------------------
 142 void AliRDHFCutsDStartoKpipi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
 143   //
 144   // Fills in vars the values of the variables
 145   //
 146   if(nvars!=fnVarsForOpt) {
 147     printf("AliRDHFCutsDStartoKpipi::GetCutsVarsForOpt: wrong number of variables\n");
 148     return;
 149   }
 150
 151
 152   AliAODRecoCascadeHF* dstarD0pi = (AliAODRecoCascadeHF*)d;
 153
 154   AliAODTrack *softPi = (AliAODTrack*)dstarD0pi->GetBachelor();
 155
 156   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)dstarD0pi->Get2Prong();
 157
 158    Int_t iter=-1;
 159   if(fVarsForOpt[0]){
 160     iter++;
 161     if(TMath::Abs(pdgdaughters[0])==211) {
 162       vars[iter]=dd->InvMassD0();
 163     } else {
 164       vars[iter]=dd->InvMassD0bar();
 165     }
 166   }
 167   if(fVarsForOpt[1]){
 168     iter++;
 169     vars[iter]=dd->GetDCA();
 170   }
 171   if(fVarsForOpt[2]){
 172     iter++;
 173     if(TMath::Abs(pdgdaughters[0])==211) {
 174       vars[iter] = dd->CosThetaStarD0();
 175     } else {
 176       vars[iter] = dd->CosThetaStarD0bar();
 177     }
 178   }
 179   if(fVarsForOpt[3]){
 180     iter++;
 181    if(TMath::Abs(pdgdaughters[0])==321) {
 182      vars[iter]=dd->PtProng(0);
 183    }
 184    else{
 185      vars[iter]=dd->PtProng(1);
 186    }
 187   }
 188   if(fVarsForOpt[4]){
 189     iter++;
 190    if(TMath::Abs(pdgdaughters[0])==211) {
 191      vars[iter]=dd->PtProng(0);
 192    }
 193    else{
 194      vars[iter]=dd->PtProng(1);
 195    }
 196   }
 197   if(fVarsForOpt[5]){
 198     iter++;
 199     if(TMath::Abs(pdgdaughters[0])==321) {
 200      vars[iter]=dd->Getd0Prong(0);
 201    }
 202    else{
 203      vars[iter]=dd->Getd0Prong(1);
 204    }
 205   }
 206   if(fVarsForOpt[6]){
 207     iter++;
 208      if(TMath::Abs(pdgdaughters[0])==211) {
 209      vars[iter]=dd->Getd0Prong(0);
 210    }
 211    else{
 212      vars[iter]=dd->Getd0Prong(1);
 213    }
 214   }
 215   if(fVarsForOpt[7]){
 216     iter++;
 217     vars[iter]= dd->Prodd0d0();
 218   }
 219   if(fVarsForOpt[8]){
 220     iter++;
 221     vars[iter]=dd->CosPointingAngle();
 222   }
 223   if(fVarsForOpt[9]){
 224     iter++;
 225     vars[iter]=dstarD0pi->InvMassDstarKpipi();
 226   }
 227   if(fVarsForOpt[10]){
 228     iter++;
 229     vars[iter]=dstarD0pi->DeltaInvMass();
 230   }
 231   if(fVarsForOpt[11]){
 232     iter++;
 233     vars[iter] = softPi->Pt();
 234   }
 235   if(fVarsForOpt[12]){
 236     iter++;
 237     vars[iter] = softPi->Pt();
 238   }
 239   if(fVarsForOpt[13]){
 240     iter++;
 241     vars[iter] =dstarD0pi->AngleD0dkpPisoft();
 242   }
 243   if(fVarsForOpt[14]){
 244     iter++;
 245     vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
 246   }
 247   if(fVarsForOpt[15]){
 248     iter++;
 249     vars[iter]=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
 250   }
 251
 252   return;
 253 }
 254 //---------------------------------------------------------------------------
 255 Int_t AliRDHFCutsDStartoKpipi::IsSelected(TObject* obj,Int_t selectionLevel) {
 256   //
 257   // Apply selection for D*.
 258   //
 259   if(!fCutsRD){
 260     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 261     return 0;
 262   }
 263
 264   AliAODRecoCascadeHF* d = (AliAODRecoCascadeHF*)obj;
 265   if(!d){
 266     cout<<"AliAODRecoCascadeHF null"<<endl;
 267     return 0;
 268   }
 269
 270   Double_t ptD=d->Pt();
 271   if(ptD<fMinPtCand) return 0;
 272   if(ptD>fMaxPtCand) return 0;
 273
 274
 275   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)d->Get2Prong();
 276   if(!dd){
 277     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 278     return 0;
 279   }
 280
 281   if(dd->HasBadDaughters()) return 0;
 282
 283   AliAODTrack *b = (AliAODTrack*)d->GetBachelor();
 284   if(fTrackCutsSoftPi && fTrackCutsSoftPi->GetRequireTPCRefit()){
 285     if(!(b->TestFilterMask(BIT(4)))) return 0;
 286   }
 287
 288   Int_t returnvalue=1;
 289   Int_t returnvaluePID=3;
 290
 291
 292   // selection on candidate
 293   if(selectionLevel==AliRDHFCuts::kAll ||
 294      selectionLevel==AliRDHFCuts::kCandidate) {
 295
 296     Double_t pt=d->Pt();
 297     Int_t ptbin=PtBin(pt);
 298
 299     // DStarMass and D0mass
 300     Double_t mDSPDG = TDatabasePDG::Instance()->GetParticle(413)->Mass();
 301     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 302     // delta mass PDG
 303     Double_t deltaPDG = mDSPDG-mD0PDG;
 304
 305     // Half width DStar mass
 306     if(TMath::Abs(mDSPDG - (d->InvMassDstarKpipi()))>fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
 307     // Half width Delta mass
 308
 309     if(TMath::Abs(deltaPDG-(d->DeltaInvMass())) > fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
 310
 311     // cut on soft pion pt
 312     if(b->Pt() < fCutsRD[GetGlobalIndex(11,ptbin)] || b->Pt() > fCutsRD[GetGlobalIndex(12,ptbin)]) return 0;
 313     // cut on the angle between D0 decay plane and soft pion
 314     if(d->AngleD0dkpPisoft() > fCutsRD[GetGlobalIndex(13,ptbin)]) return 0;
 315
 316     // select D0 that passes D* cuts
 317     returnvalue = IsD0FromDStarSelected(pt,dd,selectionLevel);
 318     if((b->Charge()==+1 && returnvalue==2) || (b->Charge()==-1 && returnvalue==1)) return 0;
 319
 320   }
 321
 322   // selection on PID
 323   if(selectionLevel==AliRDHFCuts::kAll ||
 324      selectionLevel==AliRDHFCuts::kCandidate ||
 325      selectionLevel==AliRDHFCuts::kPID) {
 326     returnvaluePID = IsSelectedPID(d);
 327   }
 328   if(returnvaluePID!=3) returnvalue =0;
 329
 330
 331   // selection on daughter tracks
 332   if(selectionLevel==AliRDHFCuts::kAll ||
 333      selectionLevel==AliRDHFCuts::kTracks) {
 334     if(!AreDaughtersSelected(dd)) return 0;
 335     if(fTrackCutsSoftPi) {
 336       AliAODVertex *vAOD = d->GetPrimaryVtx();
 337       Double_t pos[3],cov[6];
 338       vAOD->GetXYZ(pos);
 339       vAOD->GetCovarianceMatrix(cov);
 340       const AliESDVertex vESD(pos,cov,100.,100);
 341       if(!IsDaughterSelected(b,&vESD,fTrackCutsSoftPi)) return 0;
 342     }
 343   }
 344
 345   return returnvalue;
 346
 347 }
 348 //_________________________________________________________________________________________________
 349 Int_t AliRDHFCutsDStartoKpipi::IsD0FromDStarSelected(Double_t pt, TObject* obj,Int_t selectionLevel) const {
 350   //
 351   // Apply selection for D0 from D*. The selection in on D0 prongs
 352   //
 353
 354   if(!fCutsRD){
 355     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 356     return 0;
 357   }
 358
 359   AliAODRecoDecayHF2Prong* dd = (AliAODRecoDecayHF2Prong*)obj;
 360
 361   if(!dd){
 362     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 363     return 0;
 364   }
 365
 366   // selection on daughter tracks is done in IsSelected()
 367
 368   Int_t returnvalue=1;
 369
 370   // selection on candidate
 371   if(selectionLevel==AliRDHFCuts::kAll ||
 372      selectionLevel==AliRDHFCuts::kCandidate) {
 373
 374     // D0 mass
 375     Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 376     // delta mass PDG
 377
 378     Int_t ptbin=PtBin(pt);
 379
 380     Double_t mD0,mD0bar,ctsD0,ctsD0bar;
 381
 382     Int_t okD0     =0;
 383     Int_t okD0bar  =0;
 384     okD0=1; okD0bar=1;
 385
 386     if(dd->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
 387     if(dd->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || dd->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
 388
 389     if(!okD0 && !okD0bar) return 0;
 390
 391     if(TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] ||
 392        TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
 393     if(TMath::Abs(dd->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
 394        TMath::Abs(dd->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
 395     if(!okD0 && !okD0bar) return 0;
 396
 397     if(dd->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) return 0;
 398
 399     dd->InvMassD0(mD0,mD0bar);
 400     if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
 401     if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0bar = 0;
 402     if(!okD0 && !okD0bar) return 0;
 403
 404     dd->CosThetaStarD0(ctsD0,ctsD0bar);
 405     if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0;
 406     if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
 407     if(!okD0 && !okD0bar) return 0;
 408
 409     if(dd->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
 410
 411     if(dd->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
 412
 413     if(TMath::Abs(dd->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(14,ptbin)])  return 0;
 414
 415     Double_t normalDecayLengXY=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
 416     if (normalDecayLengXY < fCutsRD[GetGlobalIndex(15, ptbin)]) return 0;
 417
 418     if (okD0) returnvalue=1; //cuts passed as D0
 419     if (okD0bar) returnvalue=2; //cuts passed as D0bar
 420     if (okD0 && okD0bar) returnvalue=3; //both
 421   }
 422
 423   return returnvalue;
 424 }
 425 //----------------------------------------------------------------------------------
 426 Bool_t AliRDHFCutsDStartoKpipi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 427 {
 428   //
 429   // D* fiducial acceptance region
 430   //
 431
 432   if(pt > 5.) {
 433     // applying cut for pt > 5 GeV
 434     AliDebug(4,Form("pt of D* = %f (> 5), cutting at |y| < 0.8\n",pt));
 435     if (TMath::Abs(y) > 0.8){
 436       return kFALSE;
 437     }
 438   } else {
 439     // appliying smooth cut for pt < 5 GeV
 440     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 441     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 442     AliDebug(2,Form("pt of D* = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 443     if (y < minFiducialY || y > maxFiducialY){
 444       return kFALSE;
 445     }
 446   }
 447
 448   return kTRUE;
 449 }
 450
 451 //_______________________________________________________________________________-
 452 Int_t AliRDHFCutsDStartoKpipi::IsSelectedPID(AliAODRecoDecayHF* obj)
 453 {
 454   //
 455   // PID method, n signa approach default
 456   //
 457
 458   if(!fUsePID) return 3;
 459
 460   AliAODRecoCascadeHF* dstar = (AliAODRecoCascadeHF*)obj;
 461   if(!dstar){
 462     cout<<"AliAODRecoCascadeHF null"<<endl;
 463     return 0;
 464   }
 465   AliAODRecoDecayHF2Prong* d0 = (AliAODRecoDecayHF2Prong*)dstar->Get2Prong();
 466   if(!d0){
 467     cout<<"AliAODRecoDecayHF2Prong null"<<endl;
 468     return 0;
 469   }
 470
 471   //  here the PID
 472   AliAODTrack *pos = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(0);
 473   AliAODTrack *neg = (AliAODTrack*)dstar->Get2Prong()->GetDaughter(1);
 474
 475   if (dstar->Charge()>0){
 476     if(!SelectPID(pos,2)) return 0;//pion+
 477     if(!SelectPID(neg,3)) return 0;//kaon-
 478   }else{
 479     if(!SelectPID(pos,3)) return 0;//kaon+
 480     if(!SelectPID(neg,2)) return 0;//pion-
 481   }
 482
 483   return 3;
 484 }
 485
 486 //_______________________________________________________________________________-
 487 Int_t AliRDHFCutsDStartoKpipi::SelectPID(AliAODTrack *track, Int_t type)
 488 {
 489   //
 490   //  here the PID
 491
 492   Bool_t isParticle=kTRUE;
 493
 494   if(fPidHF->GetMatch()==1){//n-sigma
 495     Bool_t TPCon=TMath::Abs(2)>1e-4?kTRUE:kFALSE;
 496     Bool_t TOFon=TMath::Abs(3)>1e-4?kTRUE:kFALSE;
 497
 498     Bool_t isTPC=kTRUE;
 499     Bool_t isTOF=kTRUE;
 500
 501     if (TPCon){//TPC
 502       if(fPidHF->CheckStatus(track,"TPC")){
 503         if(type==2) isTPC=fPidHF->IsPionRaw(track,"TPC");
 504         if(type==3) isTPC=fPidHF->IsKaonRaw(track,"TPC");
 505       }
 506     }
 507     if (TOFon){//TOF
 508       if(fPidHF->CheckStatus(track,"TOF")){
 509         if(type==2) isTOF=fPidHF->IsPionRaw(track,"TOF");
 510         if(type==3) isTOF=fPidHF->IsKaonRaw(track,"TOF");
 511       }
 512     }
 513     isParticle = isTPC&&isTOF;
 514   }
 515
 516   if(fPidHF->GetMatch()==2){//bayesian
 517     //Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 518     Double_t prob[5]={1.,1.,1.,1.,1.};
 519
 520     //fPidHF->SetPriors(priors);
 521     fPidHF->BayesianProbability(track,prob);
 522
 523     Double_t max=0.;
 524     Int_t k=-1;
 525     for (Int_t i=0; i<5; i++) {
 526       if (prob[i]>max) {k=i; max=prob[i];}
 527     }
 528     isParticle = Bool_t(k==type);
 529   }
 530
 531   return isParticle;
 532
 533 }
 534 //__________________________________________________________________________________-
 535 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPP2010() {
 536   //
 537   //STANDARD CUTS USED FOR 2010 pp analysis
 538   //
 539   // Need to be updated for the final cut version
 540   //
 541
 542   SetName("DStartoD0piCutsStandard");
 543   SetTitle("Standard Cuts for D* analysis");
 544
 545   // PILE UP REJECTION
 546   SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
 547
 548   // EVENT CUTS
 549   SetMinVtxContr(1);
 550
 551   // CUTS ON SINGLE TRACKS
 552   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 553   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 554   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 555   esdTrackCuts->SetRequireITSRefit(kTRUE);
 556   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 557   esdTrackCuts->SetMinDCAToVertexXY(0.);
 558   esdTrackCuts->SetEtaRange(-0.8,0.8);
 559   esdTrackCuts->SetPtRange(0.3,1.e10);
 560
 561   // CUTS on SOFT PION
 562   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
 563   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
 564   esdSoftPicuts->SetRequireTPCRefit(kFALSE);
 565   esdSoftPicuts->SetRequireITSRefit(kFALSE);
 566   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 567                                           AliESDtrackCuts::kAny);
 568   esdSoftPicuts->SetPtRange(0.0,1.e10);
 569
 570   AddTrackCuts(esdTrackCuts);
 571   AddTrackCutsSoftPi(esdSoftPicuts);
 572
 573   const Int_t nptbins =13;
 574   const Double_t ptmax = 9999.;
 575   const Int_t nvars=16;
 576   Float_t ptbins[nptbins+1];
 577   ptbins[0]=0.;
 578   ptbins[1]=0.5;
 579   ptbins[2]=1.;
 580   ptbins[3]=2.;
 581   ptbins[4]=3.;
 582   ptbins[5]=4.;
 583   ptbins[6]=5.;
 584   ptbins[7]=6.;
 585   ptbins[8]=7.;
 586   ptbins[9]=8.;
 587   ptbins[10]=12.;
 588   ptbins[11]=16.;
 589   ptbins[12]=24.;
 590   ptbins[13]=ptmax;
 591
 592   SetGlobalIndex(nvars,nptbins);
 593   SetPtBins(nptbins+1,ptbins);
 594
 595   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*/
 596                                                   {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*/
 597                                                   {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 */
 598                                                   {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 */
 599                                                   {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 */
 600                                                   {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 */
 601                                                   {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 */
 602                                                   {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 */
 603                                                   {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 */
 604                                                   {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 */
 605                                                   {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 */
 606                                                   {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 */
 607                                                   {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 */
 608
 609
 610   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 611   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 612   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 613
 614   for (Int_t ibin=0;ibin<nptbins;ibin++){
 615     for (Int_t ivar = 0; ivar<nvars; ivar++){
 616       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
 617     }
 618   }
 619
 620   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 621
 622   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
 623   delete [] cutsMatrixTransposeStand;
 624   cutsMatrixTransposeStand=NULL;
 625
 626   // PID SETTINGS FOR D* analysis
 627   AliAODPidHF* pidObj=new AliAODPidHF();
 628   //pidObj->SetName("pid4DSatr");
 629   Int_t mode=1;
 630   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 631   pidObj->SetPriors(priors);
 632   pidObj->SetMatch(mode);
 633   pidObj->SetSigma(0,2); // TPC
 634   pidObj->SetSigma(3,3); // TOF
 635   pidObj->SetTPC(kTRUE);
 636   pidObj->SetTOF(kTRUE);
 637
 638   SetPidHF(pidObj);
 639   SetUsePID(kTRUE);
 640
 641   PrintAll();
 642
 643   delete pidObj;
 644   pidObj=NULL;
 645
 646   return;
 647 }
 648 //_____________________________________________________________________________-
 649 void  AliRDHFCutsDStartoKpipi::SetStandardCutsPbPb2010(){
 650   //
 651   // TEMPORARY, WORK IN PROGRESS ... BUT WORKING!
 652   //
 653   //  Lead Lead
 654   //
 655
 656   SetName("DStartoD0piCutsStandard");
 657   SetTitle("Standard Cuts for D* analysis in PbPb 2010");
 658
 659   // EVENT CUTS
 660   SetMinVtxContr(1);
 661
 662   // CUTS ON SINGLE TRACKS
 663   AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 664   esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 665   esdTrackCuts->SetRequireTPCRefit(kTRUE);
 666   esdTrackCuts->SetRequireITSRefit(kTRUE);
 667   esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 668   esdTrackCuts->SetMinDCAToVertexXY(0.);
 669   esdTrackCuts->SetEtaRange(-0.8,0.8);
 670   esdTrackCuts->SetPtRange(0.3,1.e10);
 671
 672   // CUTS on SOFT PION
 673   AliESDtrackCuts* esdSoftPicuts=new AliESDtrackCuts();
 674   esdSoftPicuts->SetRequireSigmaToVertex(kFALSE);
 675   esdSoftPicuts->SetRequireTPCRefit(kTRUE);
 676   esdSoftPicuts->SetRequireITSRefit(kTRUE);
 677   esdSoftPicuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 678                                           AliESDtrackCuts::kAny); //test d0 asimmetry
 679   esdSoftPicuts->SetPtRange(0.25,5);
 680
 681   AddTrackCuts(esdTrackCuts);
 682   AddTrackCutsSoftPi(esdSoftPicuts);
 683
 684   const Int_t nptbins =13;
 685   const Double_t ptmax = 9999.;
 686   const Int_t nvars=16;
 687   Float_t ptbins[nptbins+1];
 688   ptbins[0]=0.;
 689   ptbins[1]=0.5;
 690   ptbins[2]=1.;
 691   ptbins[3]=2.;
 692   ptbins[4]=3.;
 693   ptbins[5]=4.;
 694   ptbins[6]=5.;
 695   ptbins[7]=6.;
 696   ptbins[8]=7.;
 697   ptbins[9]=8.;
 698   ptbins[10]=12.;
 699   ptbins[11]=16.;
 700   ptbins[12]=24.;
 701   ptbins[13]=ptmax;
 702
 703   SetGlobalIndex(nvars,nptbins);
 704   SetPtBins(nptbins+1,ptbins);
 705
 706   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*/
 707                                                   {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*/
 708                                                   {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 */
 709                                                   {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 */
 710                                                   {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 */
 711                                                   {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 */
 712                                                   {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 */
 713                                                   {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 */
 714                                                   {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 */
 715                                                   {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 */
 716                                                   {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 */
 717                                                   {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 */
 718                                                   {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 */
 719
 720
 721   //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
 722   Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
 723   for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
 724
 725   for (Int_t ibin=0;ibin<nptbins;ibin++){
 726     for (Int_t ivar = 0; ivar<nvars; ivar++){
 727       cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];
 728     }
 729   }
 730
 731   SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
 732
 733   for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
 734   delete [] cutsMatrixTransposeStand;
 735   cutsMatrixTransposeStand=NULL;
 736
 737   // PID SETTINGS
 738   AliAODPidHF* pidObj=new AliAODPidHF();
 739   // pidObj->SetName("pid4DSatr");
 740   Int_t mode=1;
 741   Double_t priors[5]={0.01,0.001,0.3,0.3,0.3};
 742   pidObj->SetPriors(priors);
 743   pidObj->SetMatch(mode);
 744   pidObj->SetSigma(0,2); // TPC
 745   pidObj->SetSigma(3,3); // TOF
 746   pidObj->SetTPC(kTRUE);
 747   pidObj->SetTOF(kTRUE);
 748
 749   SetPidHF(pidObj);
 750   SetUsePID(kTRUE);
 751
 752   PrintAll();
 753
 754   delete pidObj;
 755   pidObj=NULL;
 756
 757   return;
 758
 759 }