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