PWG3/vertexingHF/AliRDHFCutsDstoKKpi.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 Ds->KKpi
  21 //
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TDatabasePDG.h>
  26 #include <Riostream.h>
  27
  28 #include "AliRDHFCutsDstoKKpi.h"
  29 #include "AliAODRecoDecayHF3Prong.h"
  30 #include "AliAODTrack.h"
  31 #include "AliESDtrack.h"
  32
  33 ClassImp(AliRDHFCutsDstoKKpi)
  34
  35 //--------------------------------------------------------------------------
  36 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const char* name) :
  37 AliRDHFCuts(name)
  38 {
  39   //
  40   // Default Constructor
  41   //
  42   Int_t nvars=14;
  43   SetNVars(nvars);
  44   TString varNames[14]={"inv. mass [GeV]",
  45                         "pTK [GeV/c]",
  46                         "pTPi [GeV/c]",
  47                         "d0K [cm]",
  48                         "d0Pi [cm]",
  49                         "dist12 [cm]",
  50                         "sigmavert [cm]",
  51                         "decLen [cm]",
  52                         "ptMax [GeV/c]",
  53                         "cosThetaPoint",
  54                         "Sum d0^2 (cm^2)",
  55                         "dca [cm]",
  56                         "inv. mass (Mphi-MKK) [GeV]",
  57                         "inv. mass (MKo*-MKpi) [GeV]"};
  58   Bool_t isUpperCut[14]={kTRUE,
  59                          kFALSE,
  60                          kFALSE,
  61                          kFALSE,
  62                          kFALSE,
  63                          kFALSE,
  64                          kTRUE,
  65                          kFALSE,
  66                          kFALSE,
  67                          kFALSE,
  68                          kFALSE,
  69                          kTRUE,
  70                          kTRUE,
  71                          kTRUE};
  72   SetVarNames(14,varNames,isUpperCut);
  73   Bool_t forOpt[14]={kFALSE,
  74                     kFALSE,
  75                     kFALSE,
  76                     kFALSE,
  77                     kFALSE,
  78                     kFALSE,
  79                     kTRUE,
  80                     kTRUE,
  81                     kTRUE,
  82                     kTRUE,
  83                     kTRUE,
  84                     kFALSE,
  85                     kTRUE,
  86                     kTRUE};
  87   SetVarsForOpt(7,forOpt);
  88   Float_t limits[2]={0,999999999.};
  89   SetPtBins(2,limits);
  90   if(fPidHF)delete fPidHF;
  91   fPidHF=new AliAODPidHF();
  92   Double_t plim[2]={0.6,0.8};
  93   Double_t nsigma[5]={2.,1.,2.,3.,0.};
  94
  95   fPidHF->SetPLimit(plim);
  96   fPidHF->SetAsym(kTRUE);
  97   fPidHF->SetSigma(nsigma);
  98   fPidHF->SetMatch(1);
  99   fPidHF->SetTPC(1);
 100   fPidHF->SetTOF(1);
 101   fPidHF->SetITS(0);
 102   fPidHF->SetTRD(0);
 103   fPidHF->SetCompat(kTRUE);
 104
 105 }
 106 //--------------------------------------------------------------------------
 107 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
 108   AliRDHFCuts(source)
 109 {
 110   //
 111   // Copy constructor
 112   //
 113
 114 }
 115 //--------------------------------------------------------------------------
 116 AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
 117 {
 118   //
 119   // assignment operator
 120   //
 121   if(&source == this) return *this;
 122
 123   AliRDHFCuts::operator=(source);
 124
 125   return *this;
 126 }
 127
 128
 129 //---------------------------------------------------------------------------
 130 void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
 131   //
 132   // Fills in vars the values of the variables
 133   //
 134
 135   if(nvars!=fnVarsForOpt) {
 136     printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
 137     return;
 138   }
 139
 140   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 141
 142   Int_t iter=-1;
 143   if(fVarsForOpt[0]){
 144     iter++;
 145     if(TMath::Abs(pdgdaughters[0]==321)){
 146       vars[iter]=dd->InvMassDsKKpi();
 147     }else{
 148       vars[iter]=dd->InvMassDspiKK();
 149     }
 150   }
 151   if(fVarsForOpt[1]){
 152     iter++;
 153     Float_t minPtDau=99999.;
 154     for(Int_t iprong=0;iprong<3;iprong++){
 155       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 156          dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
 157     }
 158     vars[iter]=minPtDau;
 159   }
 160   if(fVarsForOpt[2]){
 161     iter++;
 162     for(Int_t iprong=0;iprong<3;iprong++){
 163       if(TMath::Abs(pdgdaughters[iprong])==211) {
 164         vars[iter]=dd->PtProng(iprong);
 165       }
 166     }
 167   }
 168   if(fVarsForOpt[3]){
 169     iter++;
 170     Float_t minImpParDau=99999.;
 171     for(Int_t iprong=0;iprong<3;iprong++){
 172       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 173          dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
 174     }
 175     vars[iter]=minImpParDau;
 176   }
 177   if(fVarsForOpt[4]){
 178     iter++;
 179     for(Int_t iprong=0;iprong<3;iprong++){
 180       if(TMath::Abs(pdgdaughters[iprong])==211) {
 181         vars[iter]=dd->Getd0Prong(iprong);
 182       }
 183     }
 184   }
 185   if(fVarsForOpt[5]){
 186     iter++;
 187     Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
 188     vars[iter]=minDistPair;
 189   }
 190   if(fVarsForOpt[6]){
 191     iter++;
 192     vars[iter]=dd->GetSigmaVert();
 193   }
 194   if(fVarsForOpt[7]){
 195     iter++;
 196     vars[iter] = dd->DecayLength();
 197   }
 198   if(fVarsForOpt[8]){
 199     iter++;
 200     Float_t ptmax=0;
 201     for(Int_t i=0;i<3;i++){
 202       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
 203     }
 204     vars[iter]=ptmax;
 205   }
 206   if(fVarsForOpt[9]){
 207     iter++;
 208     vars[iter]=dd->CosPointingAngle();
 209   }
 210   if(fVarsForOpt[10]){
 211     iter++;
 212     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
 213   }
 214   if(fVarsForOpt[11]){
 215     iter++;
 216     Float_t maxDCA=0.;
 217     for(Int_t i=0;i<3;i++){
 218       if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
 219     }
 220     vars[iter]=maxDCA;
 221   }
 222   if(fVarsForOpt[12]){
 223     iter++;
 224     Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 225     if(TMath::Abs(pdgdaughters[0]==321)){
 226
 227       Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
 228        vars[iter]=TMath::Abs(phimass01-mPDGPhi);
 229        // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
 230     }else{
 231       Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
 232        vars[iter]=TMath::Abs(phimass12-mPDGPhi);
 233        // vars[iter]=dd->InvMass2Prongs(1,2,321,321);
 234     }
 235   }
 236   if(fVarsForOpt[13]){
 237     iter++;
 238     Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 239     if(TMath::Abs(pdgdaughters[0]==321)){
 240
 241       Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
 242       vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
 243       //              vars[iter]=dd->InvMass2Prongs(1,2,321,211);
 244     }else{
 245       Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
 246       vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
 247       //        vars[iter]=dd->InvMass2Prongs(0,1,211,321);
 248     }
 249   }
 250
 251
 252   return;
 253 }
 254 //---------------------------------------------------------------------------
 255 Bool_t AliRDHFCutsDstoKKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 256 {
 257   //
 258   // Checking if Ds is in fiducial acceptance region
 259   //
 260
 261   if(pt > 5.) {
 262     // applying cut for pt > 5 GeV
 263     AliDebug(2,Form("pt of Ds = %f (> 5), cutting at |y| < 0.8",pt));
 264     if (TMath::Abs(y) > 0.8) return kFALSE;
 265
 266   } else {
 267     // appliying smooth cut for pt < 5 GeV
 268     Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 269     Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 270     AliDebug(2,Form("pt of Ds = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 271     if (y < minFiducialY || y > maxFiducialY) return kFALSE;
 272   }
 273
 274   return kTRUE;
 275 }
 276
 277 //---------------------------------------------------------------------------
 278 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
 279   // PID selection
 280   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both
 281   Int_t retCode=3;
 282   Bool_t okKKpi=kTRUE;
 283   Bool_t okpiKK=kTRUE;
 284   if(!fUsePID || !rd) return retCode;
 285   if(!fPidHF){
 286     AliWarning("AliAODPidHF not created!");
 287     return retCode;
 288   }
 289   Int_t nKaons=0;
 290   Int_t nNotKaons=0;
 291   Int_t sign= rd->GetCharge();
 292   for(Int_t iDaught=0; iDaught<3; iDaught++){
 293     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
 294     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 295     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 296     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
 297
 298     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
 299     if(sign!=track->Charge()){// must be kaon
 300       if(isKaon<0) return 0;
 301     }
 302     if(isKaon>0 && isPion<0) nKaons++;
 303     if(isKaon<0) nNotKaons++;
 304     if(iDaught==0){
 305       if(isKaon<0) okKKpi=kFALSE;
 306       if(isPion<0) okpiKK=kFALSE;
 307     }
 308     else if(iDaught==2){
 309       if(isKaon<0) okpiKK=kFALSE;
 310       if(isPion<0) okKKpi=kFALSE;
 311     }
 312   }
 313
 314   if(nKaons>2)return 0;
 315   if(nNotKaons>1) return 0;
 316
 317   if(!okKKpi) retCode-=1;
 318   if(!okpiKK) retCode-=2;
 319
 320   return retCode;
 321 }
 322
 323 //---------------------------------------------------------------------------
 324 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel, AliAODEvent* aod) {
 325   //
 326   // Apply selection
 327   //
 328
 329   if(!fCutsRD){
 330     cout<<"Cut matrix not inizialized. Exit..."<<endl;
 331     return 0;
 332   }
 333   //PrintAll();
 334   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 335
 336   if(!d){
 337     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
 338     return 0;
 339   }
 340
 341   Double_t ptD=d->Pt();
 342   if(ptD<fMinPtCand) return 0;
 343   if(ptD>fMaxPtCand) return 0;
 344
 345   // selection on daughter tracks
 346   if(selectionLevel==AliRDHFCuts::kAll ||
 347      selectionLevel==AliRDHFCuts::kTracks) {
 348     if(!AreDaughtersSelected(d)) return 0;
 349   }
 350
 351
 352
 353
 354   // selection on candidate
 355   if(selectionLevel==AliRDHFCuts::kAll ||
 356      selectionLevel==AliRDHFCuts::kCandidate) {
 357     //recalculate vertex w/o daughters
 358     AliAODVertex *origownvtx=0x0;
 359     if(fRemoveDaughtersFromPrimary) {
 360       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
 361       if(!RecalcOwnPrimaryVtx(d,aod)) {
 362         CleanOwnPrimaryVtx(d,aod,origownvtx);
 363         return 0;
 364       }
 365     }
 366
 367     Int_t okDsKKpi=1;
 368     Int_t okDspiKK=1;
 369     Int_t okMassPhi=0;
 370     Int_t okMassK0star=0;
 371
 372     Double_t pt=d->Pt();
 373     Int_t ptbin=PtBin(pt);
 374     if (ptbin==-1) {
 375       CleanOwnPrimaryVtx(d,aod,origownvtx);
 376       return 0;
 377     }
 378
 379     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
 380     Double_t mDsKKpi=d->InvMassDsKKpi();
 381     Double_t mDspiKK=d->InvMassDspiKK();
 382     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
 383     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
 384     if(!okDsKKpi && !okDspiKK){
 385       CleanOwnPrimaryVtx(d,aod,origownvtx);
 386       return 0;
 387     }
 388
 389
 390
 391     // cuts on resonant decays (via Phi or K0*)
 392     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 393     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 394     if(okDsKKpi){
 395       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
 396       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
 397       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
 398       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
 399       if(!okMassPhi && !okMassK0star) okDsKKpi=0;
 400     }
 401     if(okDspiKK){
 402       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
 403       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
 404       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
 405       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
 406       if(!okMassPhi && !okMassK0star) okDspiKK=0;
 407     }
 408     if(!okDsKKpi && !okDspiKK){
 409       CleanOwnPrimaryVtx(d,aod,origownvtx);
 410       return 0;
 411     }
 412
 413     // Cuts on track pairs
 414     for(Int_t i=0;i<3;i++){
 415       if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]){
 416         CleanOwnPrimaryVtx(d,aod,origownvtx);
 417         return 0;
 418       }
 419     }
 420     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] ||
 421        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]){
 422       CleanOwnPrimaryVtx(d,aod,origownvtx);
 423       return 0;
 424     }
 425
 426
 427
 428     //single track
 429     if(TMath::Abs(d->Pt2Prong(1)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 430        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]){
 431       CleanOwnPrimaryVtx(d,aod,origownvtx);
 432       return 0;
 433     }
 434
 435     if(okDsKKpi){
 436       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 437          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
 438       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 439          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
 440     }
 441     if(okDspiKK){
 442       if(TMath::Abs(d->Pt2Prong(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]*fCutsRD[GetGlobalIndex(2,ptbin)] ||
 443          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
 444       if(TMath::Abs(d->Pt2Prong(2)) < fCutsRD[GetGlobalIndex(1,ptbin)]*fCutsRD[GetGlobalIndex(1,ptbin)] ||
 445          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
 446     }
 447     if(!okDsKKpi && !okDspiKK){
 448       CleanOwnPrimaryVtx(d,aod,origownvtx);
 449       return 0;
 450     }
 451
 452     // Cuts on candidate triplet
 453
 454     if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]){
 455       CleanOwnPrimaryVtx(d,aod,origownvtx);
 456       return 0;
 457     }
 458
 459     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]){
 460       CleanOwnPrimaryVtx(d,aod,origownvtx);
 461       return 0;
 462     }
 463
 464     if(d->Pt2Prong(0)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 465        d->Pt2Prong(1)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)] &&
 466        d->Pt2Prong(2)<fCutsRD[GetGlobalIndex(8,ptbin)]*fCutsRD[GetGlobalIndex(8,ptbin)]) {
 467       CleanOwnPrimaryVtx(d,aod,origownvtx);
 468       return 0;
 469     }
 470
 471     if(d->DecayLength2()<fCutsRD[GetGlobalIndex(7,ptbin)]*fCutsRD[GetGlobalIndex(7,ptbin)]){
 472       CleanOwnPrimaryVtx(d,aod,origownvtx);
 473       return 0;
 474     }
 475
 476
 477     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 478     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)]){
 479       CleanOwnPrimaryVtx(d,aod,origownvtx);
 480       return 0;
 481     }
 482
 483      // unset recalculated primary vertex when not needed any more
 484     CleanOwnPrimaryVtx(d,aod,origownvtx);
 485
 486
 487     // PID selection
 488     Int_t returnvaluePID=3;
 489     if(selectionLevel==AliRDHFCuts::kAll ||
 490        selectionLevel==AliRDHFCuts::kCandidate ||
 491        selectionLevel==AliRDHFCuts::kPID) {
 492       returnvaluePID = IsSelectedPID(d);
 493       fIsSelectedPID=returnvaluePID;
 494     }
 495     if(returnvaluePID==0)return 0;
 496
 497     Bool_t okPidDsKKpi=returnvaluePID&1;
 498     Bool_t okPidDspiKK=returnvaluePID&2;
 499     if((okPidDsKKpi && okDsKKpi)||(okPidDspiKK && okDspiKK)){
 500       Int_t returnvalue=0;
 501       if(okDsKKpi) returnvalue+=1;
 502       if(okDspiKK) returnvalue+=2;
 503       if(okMassPhi) returnvalue+=4;
 504       if(okMassK0star) returnvalue+=8;
 505       return returnvalue;
 506     }else{
 507       return 0;
 508     }
 509   }
 510   return 15;
 511
 512 }
 513 //---------------------------------------------------------------------------