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