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 }
  89 //--------------------------------------------------------------------------
  90 AliRDHFCutsDstoKKpi::AliRDHFCutsDstoKKpi(const AliRDHFCutsDstoKKpi &source) :
  91   AliRDHFCuts(source)
  92 {
  93   //
  94   // Copy constructor
  95   //
  96
  97 }
  98 //--------------------------------------------------------------------------
  99 AliRDHFCutsDstoKKpi &AliRDHFCutsDstoKKpi::operator=(const AliRDHFCutsDstoKKpi &source)
 100 {
 101   //
 102   // assignment operator
 103   //
 104   if(&source == this) return *this;
 105
 106   AliRDHFCuts::operator=(source);
 107
 108   return *this;
 109 }
 110
 111
 112 //---------------------------------------------------------------------------
 113 void AliRDHFCutsDstoKKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
 114   //
 115   // Fills in vars the values of the variables
 116   //
 117
 118   if(nvars!=fnVarsForOpt) {
 119     printf("AliRDHFCutsDstoKKpi::GetCutsVarsForOpt: wrong number of variables\n");
 120     return;
 121   }
 122
 123   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 124
 125   Int_t iter=-1;
 126   if(fVarsForOpt[0]){
 127     iter++;
 128     if(TMath::Abs(pdgdaughters[0]==321)){
 129       vars[iter]=dd->InvMassDsKKpi();
 130     }else{
 131       vars[iter]=dd->InvMassDspiKK();
 132     }
 133   }
 134   if(fVarsForOpt[1]){
 135     iter++;
 136     Float_t minPtDau=99999.;
 137     for(Int_t iprong=0;iprong<3;iprong++){
 138       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 139          dd->PtProng(iprong)<minPtDau) minPtDau=dd->PtProng(iprong);
 140     }
 141     vars[iter]=minPtDau;
 142   }
 143   if(fVarsForOpt[2]){
 144     iter++;
 145     for(Int_t iprong=0;iprong<3;iprong++){
 146       if(TMath::Abs(pdgdaughters[iprong])==211) {
 147         vars[iter]=dd->PtProng(iprong);
 148       }
 149     }
 150   }
 151   if(fVarsForOpt[3]){
 152     iter++;
 153     Float_t minImpParDau=99999.;
 154     for(Int_t iprong=0;iprong<3;iprong++){
 155       if(TMath::Abs(pdgdaughters[iprong])==321 &&
 156          dd->Getd0Prong(iprong)<minImpParDau) minImpParDau=dd->Getd0Prong(iprong);
 157     }
 158     vars[iter]=minImpParDau;
 159   }
 160   if(fVarsForOpt[4]){
 161     iter++;
 162     for(Int_t iprong=0;iprong<3;iprong++){
 163       if(TMath::Abs(pdgdaughters[iprong])==211) {
 164         vars[iter]=dd->Getd0Prong(iprong);
 165       }
 166     }
 167   }
 168   if(fVarsForOpt[5]){
 169     iter++;
 170     Float_t minDistPair=TMath::Min(dd->GetDist12toPrim(),dd->GetDist23toPrim());
 171     vars[iter]=minDistPair;
 172   }
 173   if(fVarsForOpt[6]){
 174     iter++;
 175     vars[iter]=dd->GetSigmaVert();
 176   }
 177   if(fVarsForOpt[7]){
 178     iter++;
 179     vars[iter] = dd->DecayLength();
 180   }
 181   if(fVarsForOpt[8]){
 182     iter++;
 183     Float_t ptmax=0;
 184     for(Int_t i=0;i<3;i++){
 185       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
 186     }
 187     vars[iter]=ptmax;
 188   }
 189   if(fVarsForOpt[9]){
 190     iter++;
 191     vars[iter]=dd->CosPointingAngle();
 192   }
 193   if(fVarsForOpt[10]){
 194     iter++;
 195     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
 196   }
 197   if(fVarsForOpt[11]){
 198     iter++;
 199     Float_t maxDCA=0.;
 200     for(Int_t i=0;i<3;i++){
 201       if(d->GetDCA(i)>maxDCA) maxDCA=d->GetDCA(i);
 202     }
 203     vars[iter]=maxDCA;
 204   }
 205   if(fVarsForOpt[12]){
 206     iter++;
 207     Double_t mPDGPhi = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 208     if(TMath::Abs(pdgdaughters[0]==321)){
 209
 210       Double_t phimass01=d->InvMass2Prongs(0,1,321,321);
 211        vars[iter]=TMath::Abs(phimass01-mPDGPhi);
 212        // vars[iter]=dd->InvMass2Prongs(0,1,321,321);
 213     }else{
 214       Double_t phimass12=d->InvMass2Prongs(1,2,321,321);
 215        vars[iter]=TMath::Abs(phimass12-mPDGPhi);
 216        // vars[iter]=dd->InvMass2Prongs(1,2,321,321);
 217     }
 218   }
 219   if(fVarsForOpt[13]){
 220     iter++;
 221     Double_t mPDGK0star = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 222     if(TMath::Abs(pdgdaughters[0]==321)){
 223
 224       Double_t mass12kpi=d->InvMass2Prongs(1,2,321,211);
 225       vars[iter]=TMath::Abs(mass12kpi-mPDGK0star);
 226       //              vars[iter]=dd->InvMass2Prongs(1,2,321,211);
 227     }else{
 228       Double_t mass01pik=d->InvMass2Prongs(0,1,211,321);
 229       vars[iter]=TMath::Abs(mass01pik-mPDGK0star);
 230       //        vars[iter]=dd->InvMass2Prongs(0,1,211,321);
 231     }
 232   }
 233
 234
 235   return;
 236 }
 237 //---------------------------------------------------------------------------
 238 Int_t AliRDHFCutsDstoKKpi::IsSelectedPID(AliAODRecoDecayHF *rd) {
 239   // PID selection
 240   // return values: 0->NOT OK, 1->OK as KKpi, 2->OK as piKK, 3->OK as both
 241   Int_t retCode=3;
 242   Bool_t okKKpi=kTRUE;
 243   Bool_t okpiKK=kTRUE;
 244   if(!fUsePID || !rd) return retCode;
 245   if(!fPidHF){
 246     AliWarning("AliAODPidHF not created!");
 247     return retCode;
 248   }
 249   Int_t nKaons=0;
 250   Int_t nNotKaons=0;
 251   Int_t sign= rd->GetCharge();
 252   for(Int_t iDaught=0; iDaught<3; iDaught++){
 253     AliAODTrack *track=(AliAODTrack*)rd->GetDaughter(iDaught);
 254     Int_t isPion=fPidHF->MakeRawPid(track,AliPID::kPion);
 255     Int_t isKaon=fPidHF->MakeRawPid(track,AliPID::kKaon);
 256     Int_t isProton=fPidHF->MakeRawPid(track,AliPID::kProton);
 257
 258     if(isProton>0 &&  isKaon<0  && isPion<0) return 0;
 259     if(sign!=track->Charge()){// must be kaon
 260       if(isKaon<0) return 0;
 261     }
 262     if(isKaon>0 && isPion<0) nKaons++;
 263     if(isKaon<0) nNotKaons++;
 264     if(iDaught==0){
 265       if(isKaon<0) okKKpi=kFALSE;
 266       if(isPion<0) okpiKK=kFALSE;
 267     }
 268     else if(iDaught==2){
 269       if(isKaon<0) okpiKK=kFALSE;
 270       if(isPion<0) okKKpi=kFALSE;
 271     }
 272   }
 273
 274   if(nKaons>2)return 0;
 275   if(nNotKaons>1) return 0;
 276
 277   if(!okKKpi) retCode-=1;
 278   if(!okpiKK) retCode-=2;
 279
 280   return retCode;
 281 }
 282
 283 //---------------------------------------------------------------------------
 284 Int_t AliRDHFCutsDstoKKpi::IsSelected(TObject* obj,Int_t selectionLevel) {
 285   //
 286   // Apply selection
 287   //
 288
 289   if(!fCutsRD){
 290     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 291     return 0;
 292   }
 293   //PrintAll();
 294   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 295
 296   if(!d){
 297     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
 298     return 0;
 299   }
 300
 301
 302   // selection on daughter tracks
 303   if(selectionLevel==AliRDHFCuts::kAll ||
 304      selectionLevel==AliRDHFCuts::kTracks) {
 305     if(!AreDaughtersSelected(d)) return 0;
 306   }
 307
 308
 309
 310   // PID selection
 311   Int_t returnvaluePID=3;
 312   if(selectionLevel==AliRDHFCuts::kAll ||
 313      selectionLevel==AliRDHFCuts::kCandidate ||
 314      selectionLevel==AliRDHFCuts::kPID) {
 315     returnvaluePID = IsSelectedPID(d);
 316   }
 317   if(returnvaluePID==0)return 0;
 318   Bool_t okPidDsKKpi=returnvaluePID&1;
 319   Bool_t okPidDspiKK=returnvaluePID&2;
 320
 321   Int_t returnvalue=1;
 322   // selection on candidate
 323   if(selectionLevel==AliRDHFCuts::kAll ||
 324      selectionLevel==AliRDHFCuts::kCandidate) {
 325
 326
 327     Int_t okDsKKpi=1;
 328     Int_t okDspiKK=1;
 329     Int_t okMassPhi=0;
 330     Int_t okMassK0star=0;
 331
 332     Double_t pt=d->Pt();
 333     Int_t ptbin=PtBin(pt);
 334
 335     Double_t mDsPDG = TDatabasePDG::Instance()->GetParticle(431)->Mass();
 336     Double_t mDsKKpi=d->InvMassDsKKpi();
 337     Double_t mDspiKK=d->InvMassDspiKK();
 338     if(TMath::Abs(mDsKKpi-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDsKKpi = 0;
 339     if(TMath::Abs(mDspiKK-mDsPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okDspiKK = 0;
 340     if(!okDsKKpi && !okDspiKK) return 0;
 341     if(okPidDsKKpi && !okDsKKpi)  return 0;
 342     if(okPidDspiKK && !okDspiKK) return 0;
 343
 344     //single track
 345     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] ||
 346        TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;
 347     if(okDsKKpi){
 348       if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(1,ptbin)] ||
 349          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDsKKpi=0;
 350       if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)] ||
 351          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDsKKpi=0;
 352     }
 353     if(okDspiKK){
 354       if(TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)] ||
 355          TMath::Abs(d->Getd0Prong(0))<fCutsRD[GetGlobalIndex(4,ptbin)]) okDspiKK=0;
 356       if(TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(1,ptbin)] ||
 357          TMath::Abs(d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(3,ptbin)]) okDspiKK=0;
 358     }
 359     if(!okDsKKpi && !okDspiKK) return 0;
 360
 361     // cuts on resonant decays (via Phi or K0*)
 362     Double_t mPhiPDG = TDatabasePDG::Instance()->GetParticle(333)->Mass();
 363     Double_t mK0starPDG = TDatabasePDG::Instance()->GetParticle(313)->Mass();
 364     if(okDsKKpi){
 365       Double_t mass01phi=d->InvMass2Prongs(0,1,321,321);
 366       Double_t mass12K0s=d->InvMass2Prongs(1,2,321,211);
 367       if(TMath::Abs(mass01phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
 368       if(TMath::Abs(mass12K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
 369       if(!okMassPhi && !okMassK0star) okDsKKpi=0;
 370     }
 371     if(okDspiKK){
 372       Double_t mass01K0s=d->InvMass2Prongs(0,1,211,321);
 373       Double_t mass12phi=d->InvMass2Prongs(1,2,321,321);
 374       if(TMath::Abs(mass01K0s-mK0starPDG)<fCutsRD[GetGlobalIndex(13,ptbin)]) okMassK0star = 1;
 375       if(TMath::Abs(mass12phi-mPhiPDG)<fCutsRD[GetGlobalIndex(12,ptbin)]) okMassPhi=1;
 376       if(!okMassPhi && !okMassK0star) okDspiKK=0;
 377     }
 378     if(!okDsKKpi && !okDspiKK) return 0;
 379
 380     // Cuts on track pairs
 381     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)])  return 0;
 382     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)] ||
 383        d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;
 384
 385
 386     // Cuts on candidate triplet
 387     if(d->GetSigmaVert()>fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;
 388     if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
 389     if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] &&
 390        TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] &&
 391        TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
 392     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)])return 0;
 393     Double_t sum2=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 394     if(sum2<fCutsRD[GetGlobalIndex(10,ptbin)])return 0;
 395
 396     returnvalue=0;
 397     if(okDsKKpi) returnvalue+=1;
 398     if(okDspiKK) returnvalue+=2;
 399     if(okMassPhi) returnvalue+=4;
 400     if(okMassK0star) returnvalue+=8;
 401
 402   }
 403   return returnvalue;
 404
 405 }
 406 //---------------------------------------------------------------------------