PWGHF/vertexingHF/AliRDHFCutsLctopKpi.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 Lc->pKpi
  21 //
  22 // Author: A.Dainese, andrea.dainese@pd.infn.it
  23 /////////////////////////////////////////////////////////////
  24
  25 #include <TDatabasePDG.h>
  26 #include <Riostream.h>
  27
  28 #include "AliRDHFCutsLctopKpi.h"
  29 #include "AliAODRecoDecayHF3Prong.h"
  30 #include "AliRDHFCuts.h"
  31 #include "AliAODTrack.h"
  32 #include "AliESDtrack.h"
  33 #include "AliKFParticle.h"
  34 #include "AliESDVertex.h"
  35
  36 ClassImp(AliRDHFCutsLctopKpi)
  37
  38 //--------------------------------------------------------------------------
  39 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const char* name) :
  40 AliRDHFCuts(name),
  41 fPidObjprot(0),
  42 fPidObjpion(0),
  43 fUseImpParProdCorrCut(kFALSE),
  44 fPIDStrategy(kNSigma),
  45 fCutsStrategy(kStandard)
  46 {
  47   //
  48   // Default Constructor
  49   //
  50   Int_t nvars=13;
  51   SetNVars(nvars);
  52   TString varNames[13]={"inv. mass [GeV]",
  53                         "pTK [GeV/c]",
  54                         "pTP [GeV/c]",
  55                         "d0K [cm]   lower limit!",
  56                         "d0Pi [cm]  lower limit!",
  57                         "dist12 (cm)",
  58                         "sigmavert (cm)",
  59                         "dist prim-sec (cm)",
  60                         "pM=Max{pT1,pT2,pT3} (GeV/c)",
  61                         "cosThetaPoint",
  62                         "Sum d0^2 (cm^2)",
  63                         "dca cut (cm)",
  64                         "cut on pTpion [GeV/c]"};
  65   Bool_t isUpperCut[13]={kTRUE,
  66                          kFALSE,
  67                          kFALSE,
  68                          kFALSE,
  69                          kFALSE,
  70                          kFALSE,
  71                          kTRUE,
  72                          kFALSE,
  73                          kFALSE,
  74                          kFALSE,
  75                          kFALSE,
  76                          kTRUE,
  77                          kFALSE
  78                          };
  79   SetVarNames(nvars,varNames,isUpperCut);
  80   Bool_t forOpt[13]={kFALSE,
  81                      kTRUE,
  82                      kTRUE,
  83                      kFALSE,
  84                      kFALSE,
  85                      kFALSE,
  86                      kFALSE,
  87                      kTRUE,
  88                      kFALSE,
  89                      kFALSE,
  90                      kFALSE,
  91                      kFALSE,
  92                      kTRUE};
  93   SetVarsForOpt(4,forOpt);
  94   Float_t limits[2]={0,999999999.};
  95   SetPtBins(2,limits);
  96   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies)
  97       fPIDThreshold[ispecies]=0.;
  98 }
  99 //--------------------------------------------------------------------------
 100 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const AliRDHFCutsLctopKpi &source) :
 101   AliRDHFCuts(source),
 102   fPidObjprot(0x0),
 103   fPidObjpion(0x0),
 104   fUseImpParProdCorrCut(source.fUseImpParProdCorrCut),
 105   fPIDStrategy(source.fPIDStrategy),
 106   fCutsStrategy(source.fCutsStrategy)
 107 {
 108   //
 109   // Copy constructor
 110   //
 111   if (source.fPidObjprot) fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
 112   else fPidObjprot = new AliAODPidHF();
 113   if (source.fPidObjpion) fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
 114   else fPidObjpion = new AliAODPidHF();
 115   memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
 116 }
 117 //--------------------------------------------------------------------------
 118 AliRDHFCutsLctopKpi &AliRDHFCutsLctopKpi::operator=(const AliRDHFCutsLctopKpi &source)
 119 {
 120   //
 121   // assignment operator
 122   //
 123   if(this != &source) {
 124
 125     AliRDHFCuts::operator=(source);
 126     delete fPidObjprot;
 127     fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
 128     delete fPidObjpion;
 129     fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
 130     fPIDStrategy=source.fPIDStrategy;
 131     fCutsStrategy=source.fCutsStrategy;
 132     memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
 133   }
 134
 135   return *this;
 136 }
 137 //---------------------------------------------------------------------------
 138 AliRDHFCutsLctopKpi::~AliRDHFCutsLctopKpi() {
 139  //
 140  //  // Default Destructor
 141  //
 142  if(fPidObjpion){
 143   delete fPidObjpion;
 144   fPidObjpion=0;
 145  }
 146  if(fPidObjprot){
 147   delete fPidObjprot;
 148   fPidObjprot=0;
 149  }
 150
 151 }
 152
 153 //---------------------------------------------------------------------------
 154 void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters, AliAODEvent *aod) {
 155   //
 156   // Fills in vars the values of the variables
 157   //
 158
 159   if(nvars!=fnVarsForOpt) {
 160     printf("AliRDHFCutsLctopKpi::GetCutsVarsForOpt: wrong number of variables\n");
 161     return;
 162   }
 163
 164   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 165
 166     Int_t iter=-1;
 167   if(fVarsForOpt[0]){
 168     iter++;
 169     vars[iter]=dd->InvMassLcpKpi();
 170   }
 171   if(fVarsForOpt[1]){
 172     iter++;
 173     for(Int_t iprong=0;iprong<3;iprong++){
 174       if(TMath::Abs(pdgdaughters[iprong])==321) {
 175         vars[iter]=dd->PtProng(iprong);
 176       }
 177     }
 178   }
 179   if(fVarsForOpt[2]){
 180     iter++;
 181     for(Int_t iprong=0;iprong<3;iprong++){
 182       if(TMath::Abs(pdgdaughters[iprong])==2212) {
 183         vars[iter]=dd->PtProng(iprong);
 184       }
 185     }
 186   }
 187   if(fVarsForOpt[3]){
 188     iter++;
 189     for(Int_t iprong=0;iprong<3;iprong++){
 190       if(TMath::Abs(pdgdaughters[iprong])==2212) {
 191         vars[iter]=dd->Getd0Prong(iprong);
 192       }
 193     }
 194   }
 195   if(fVarsForOpt[4]){
 196     iter++;
 197     for(Int_t iprong=0;iprong<3;iprong++){
 198       if(TMath::Abs(pdgdaughters[iprong])==211) {
 199         vars[iter]=dd->Getd0Prong(iprong);
 200       }
 201     }
 202   }
 203   if(fVarsForOpt[5]){
 204     iter++;
 205     vars[iter]=dd->GetDist12toPrim();
 206   }
 207   if(fVarsForOpt[6]){
 208     iter++;
 209     vars[iter]=dd->GetSigmaVert(aod);
 210   }
 211   if(fVarsForOpt[7]){
 212     iter++;
 213     vars[iter] = dd->DecayLength();
 214   }
 215   if(fVarsForOpt[8]){
 216     iter++;
 217     Float_t ptmax=0;
 218     for(Int_t i=0;i<3;i++){
 219       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
 220     }
 221     vars[iter]=ptmax;
 222   }
 223   if(fVarsForOpt[9]){
 224     iter++;
 225     vars[iter]=dd->CosPointingAngle();
 226   }
 227   if(fVarsForOpt[10]){
 228     iter++;
 229     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
 230   }
 231   if(fVarsForOpt[11]){
 232     iter++;
 233     vars[iter]=dd->GetDCA();
 234   }
 235   if(fVarsForOpt[12]){
 236     iter++;
 237     for(Int_t iprong=0;iprong<3;iprong++){
 238       if(TMath::Abs(pdgdaughters[iprong])==211) {
 239         vars[iter]=dd->PtProng(iprong);
 240       }
 241     }
 242   }
 243
 244   return;
 245 }
 246 //---------------------------------------------------------------------------
 247 Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent *aod) {
 248   //
 249   // Apply selection
 250   //
 251
 252   if(!fCutsRD){
 253     cout<<"Cut matrice not inizialized. Exit..."<<endl;
 254     return 0;
 255   }
 256   //PrintAll();
 257   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 258
 259   if(!d){
 260     cout<<"AliAODRecoDecayHF3Prong null"<<endl;
 261     return 0;
 262   }
 263
 264
 265   if(fKeepSignalMC) if(IsSignalMC(d,aod,4122)) return 3;
 266
 267   Int_t returnvalue=3;
 268   Int_t returnvaluePID=3;
 269
 270   if(d->Pt()<fMinPtCand) return 0;
 271   if(d->Pt()>fMaxPtCand) return 0;
 272
 273   if(d->HasBadDaughters()) return 0;
 274
 275   if(selectionLevel==AliRDHFCuts::kAll ||
 276      selectionLevel==AliRDHFCuts::kCandidate||
 277      selectionLevel==AliRDHFCuts::kPID) {
 278      switch (fPIDStrategy) {
 279       case kNSigma:
 280        returnvaluePID = IsSelectedPID(d);
 281
 282       break;
 283       case kCombined:
 284        returnvaluePID = IsSelectedCombinedPID(d);
 285       break;
 286      }
 287      fIsSelectedPID=returnvaluePID;
 288   }
 289   //  if(fUsePID || selectionLevel==AliRDHFCuts::kPID) returnvaluePID = IsSelectedCombinedPID(d);   // to test!!
 290   if(returnvaluePID==0) return 0;
 291
 292   // selection on candidate
 293   if(selectionLevel==AliRDHFCuts::kAll ||
 294      selectionLevel==AliRDHFCuts::kCandidate) {
 295
 296     Double_t pt=d->Pt();
 297
 298     Int_t ptbin=PtBin(pt);
 299
 300     Double_t mLcpKpi,mLcpiKp;
 301     Int_t okLcpKpi=1,okLcpiKp=1;
 302
 303     Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 304
 305     mLcpKpi=d->InvMassLcpKpi();
 306     mLcpiKp=d->InvMassLcpiKp();
 307
 308     if(TMath::Abs(mLcpKpi-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpKpi = 0;
 309     if(TMath::Abs(mLcpiKp-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpiKp = 0;
 310     if(!okLcpKpi && !okLcpiKp) return 0;
 311
 312   switch (fCutsStrategy) {
 313
 314     case kStandard:
 315     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;//Kaon
 316     if(d->Pt()>=3. && d->PProng(1)<0.55) return 0;
 317     if((TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(12,ptbin)])) okLcpKpi=0;
 318     if((TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(12,ptbin)]))okLcpiKp=0;
 319     if(!okLcpKpi && !okLcpiKp) return 0;
 320     //2track cuts
 321     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;
 322     if(d->GetDist12toPrim()>0.5) return 0;
 323     if(d->GetDist23toPrim()>0.5) return 0;
 324     if(fUseImpParProdCorrCut){
 325       if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.) return 0;
 326     }
 327     //sec vert
 328     if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
 329     if(d->DecayLength()>0.5) return 0;
 330
 331   //  Double_t sumd0s=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 332   //  if(sumd0s<fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
 333     if((d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2))<fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
 334
 335     if(TMath::Abs(d->PtProng(0))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(1))<fCutsRD[GetGlobalIndex(8,ptbin)] && TMath::Abs(d->PtProng(2))<fCutsRD[GetGlobalIndex(8,ptbin)]) return 0;
 336     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
 337     if(d->GetSigmaVert(aod)>fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;
 338
 339     //DCA
 340     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) return 0;
 341
 342     break;
 343
 344    case kKF:
 345     Int_t pdgs[3]={0,321,0};
 346     Bool_t constraint=kFALSE;
 347     if(fCutsRD[GetGlobalIndex(1,ptbin)]>0.) constraint=kTRUE;
 348     Double_t field=aod->GetMagneticField();
 349     if (returnvaluePID==1 || returnvaluePID==3){
 350
 351       pdgs[0]=2122;pdgs[2]=211;
 352       AliKFParticle *lc1=ReconstructKF(d,pdgs,field,constraint);
 353       if(!lc1){
 354         okLcpKpi=0;
 355       }else{
 356         if(lc1->GetChi2()/lc1->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)]) okLcpKpi=0;
 357       }
 358     } else if(returnvaluePID>=2){
 359
 360       pdgs[0]=211;pdgs[2]=2212;
 361       AliKFParticle *lc2=ReconstructKF(d,pdgs,field,constraint);
 362       if(!lc2){
 363         okLcpiKp=0;
 364       }else{
 365         if(lc2->GetChi2()/lc2->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)])okLcpiKp=0;
 366       }
 367     }
 368     break;
 369
 370    }
 371
 372     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 373     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 374     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 375
 376   }
 377
 378   // selection on daughter tracks
 379   if(selectionLevel==AliRDHFCuts::kAll ||
 380      selectionLevel==AliRDHFCuts::kTracks) {
 381     if(!AreDaughtersSelected(d)) return 0;
 382   }
 383
 384
 385   Int_t returnvalueTot=CombinePIDCuts(returnvalue,returnvaluePID);
 386   return returnvalueTot;
 387 }
 388 //---------------------------------------------------------------------------
 389 Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
 390
 391
 392     if(!fUsePID || !obj) return 3;
 393     Int_t okLcpKpi=0,okLcpiKp=0;
 394     Int_t returnvalue=0;
 395     Bool_t isPeriodd=fPidHF->GetOnePad();
 396     Bool_t isMC=fPidHF->GetMC();
 397     Bool_t ispion0=kTRUE,ispion2=kTRUE;
 398     Bool_t isproton0=kFALSE,isproton2=kFALSE;
 399     Bool_t iskaon1=kFALSE;
 400     if(isPeriodd) {
 401      fPidObjprot->SetOnePad(kTRUE);
 402      fPidObjpion->SetOnePad(kTRUE);
 403     }
 404     if(isMC) {
 405      fPidObjprot->SetMC(kTRUE);
 406      fPidObjpion->SetMC(kTRUE);
 407     }
 408
 409     for(Int_t i=0;i<3;i++){
 410      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
 411      if(!track) return 0;
 412      // identify kaon
 413      if(track->P()<0.55){
 414       fPidHF->SetTOF(kFALSE);
 415       fPidHF->SetTOFdecide(kFALSE);
 416      }
 417      if(i==1) {
 418       Int_t isKaon=fPidHF->MakeRawPid(track,3);
 419       if(isKaon>=1) {
 420        iskaon1=kTRUE;
 421       if(fPidHF->MakeRawPid(track,2)>=1) iskaon1=kFALSE;
 422       }
 423       if(!iskaon1) return 0;
 424
 425      }else{
 426      //pion or proton
 427     if(track->P()<1.){
 428       fPidObjprot->SetTOF(kFALSE);
 429       fPidObjprot->SetTOFdecide(kFALSE);
 430      }
 431
 432      Int_t isProton=fPidObjprot->MakeRawPid(track,4);
 433      if(isProton>=1){
 434       if(fPidHF->MakeRawPid(track,2)>=1) isProton=-1;
 435       if(fPidHF->MakeRawPid(track,3)>=1) isProton=-1;
 436      }
 437
 438      Int_t isPion=fPidObjpion->MakeRawPid(track,2);
 439      if(fPidHF->MakeRawPid(track,3)>=1) isPion=-1;
 440      if(fPidObjprot->MakeRawPid(track,4)>=1) isPion=-1;
 441      if(track->P()<1.){
 442       fPidObjprot->SetTOF(kTRUE);
 443       fPidObjprot->SetTOFdecide(kTRUE);
 444      }
 445      if(track->P()<0.55){
 446       fPidHF->SetTOF(kTRUE);
 447       fPidHF->SetTOFdecide(kTRUE);
 448      }
 449
 450      if(i==0) {
 451       if(isPion<0) ispion0=kFALSE;
 452       if(isProton>=1) isproton0=kTRUE;
 453
 454      }
 455       if(!ispion0 && !isproton0) return 0;
 456      if(i==2) {
 457       if(isPion<0) ispion2=kFALSE;
 458       if(isProton>=1) isproton2=kTRUE;
 459      }
 460
 461     }
 462    }
 463
 464     if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
 465     if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
 466     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 467     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 468     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 469
 470  return returnvalue;
 471 }
 472 //---------------------------------------------------------------------------
 473 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
 474
 475   //  Printf(" -------- IsSelectedCombinedPID --------------");
 476
 477
 478     if(!fUsePID || !obj) {return 3;}
 479     Int_t okLcpKpi=0,okLcpiKp=0;
 480     Int_t returnvalue=0;
 481     Bool_t isPeriodd=fPidHF->GetOnePad();
 482     Bool_t isMC=fPidHF->GetMC();
 483
 484     if(isPeriodd) {
 485             fPidObjprot->SetOnePad(kTRUE);
 486             fPidObjpion->SetOnePad(kTRUE);
 487     }
 488     if(isMC) {
 489             fPidObjprot->SetMC(kTRUE);
 490             fPidObjpion->SetMC(kTRUE);
 491     }
 492
 493     AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
 494     AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
 495     AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
 496     if (!track0 || !track1 || !track2) return 0;
 497     Double_t prob0[AliPID::kSPECIES];
 498     Double_t prob1[AliPID::kSPECIES];
 499     Double_t prob2[AliPID::kSPECIES];
 500     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 501     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
 502     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 503
 504    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 505     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
 506    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 507
 508     if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 509     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
 510    if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 511
 512     if(fPIDThreshold[AliPID::kPion]>0. && fPIDThreshold[AliPID::kKaon]>0. && fPIDThreshold[AliPID::kProton]>0.){
 513     okLcpiKp=  (prob0[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ])
 514              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
 515              &&(prob2[AliPID::kProton]>fPIDThreshold[AliPID::kProton]);
 516     okLcpKpi=  (prob0[AliPID::kProton]>fPIDThreshold[AliPID::kProton])
 517              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
 518              &&(prob2[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ]);
 519    }else{
 520                     //pion or proton
 521
 522
 523     if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
 524     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;
 525     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1;
 526             }
 527            }
 528
 529     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 530     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 531     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 532
 533     return returnvalue;
 534 }
 535 //-----------------------
 536 Int_t AliRDHFCutsLctopKpi::CombinePIDCuts(Int_t returnvalue, Int_t returnvaluePID) const {
 537
 538  Int_t returnvalueTot=0;
 539  Int_t okLcpKpi=0,okLcpiKp=0;
 540  if(returnvaluePID==1){
 541    if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
 542  }
 543  if(returnvaluePID==2){
 544    if(returnvalue>=2) okLcpiKp=1;
 545  }
 546  if(returnvaluePID==3 && returnvalue>0){
 547   if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
 548   if(returnvalue>=2) okLcpiKp=1;
 549  }
 550
 551  if(okLcpKpi) returnvalueTot=1; //cuts passed as Lc->pKpi
 552  if(okLcpiKp) returnvalueTot=2; //cuts passed as Lc->piKp
 553  if(okLcpKpi && okLcpiKp) returnvalueTot=3; //cuts passed as both pKpi and piKp
 554  return returnvalueTot;
 555 }
 556 //----------------------------------
 557 void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
 558
 559  SetName("LctopKpiProdCuts");
 560  SetTitle("Production cuts for Lc analysis");
 561
 562  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 563  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 564  esdTrackCuts->SetRequireTPCRefit(kTRUE);
 565  esdTrackCuts->SetMinNClustersTPC(70);
 566  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 567                                           AliESDtrackCuts::kAny);
 568  esdTrackCuts->SetRequireITSRefit(kTRUE);
 569  esdTrackCuts->SetMinNClustersITS(4);
 570  esdTrackCuts->SetMinDCAToVertexXY(0.);
 571  esdTrackCuts->SetEtaRange(-0.8,0.8);
 572  esdTrackCuts->SetPtRange(0.3,1.e10);
 573  AddTrackCuts(esdTrackCuts);
 574
 575  const Int_t nptbins=4;
 576  const Int_t nvars=13;
 577  Float_t* ptbins;
 578  ptbins=new Float_t[nptbins+1];
 579
 580  ptbins[0]=0.;
 581  ptbins[1]=2.;
 582  ptbins[2]=3.;
 583  ptbins[3]=4.;
 584  ptbins[4]=9999.;
 585
 586  SetGlobalIndex(nvars,nptbins);
 587  SetPtBins(nptbins+1,ptbins);
 588
 589  Float_t** prodcutsval;
 590  prodcutsval=new Float_t*[nvars];
 591  for(Int_t iv=0;iv<nvars;iv++){
 592   prodcutsval[iv]=new Float_t[nptbins];
 593  }
 594
 595  for(Int_t ipt=0;ipt<nptbins;ipt++){
 596   prodcutsval[0][ipt]=0.18;
 597   prodcutsval[1][ipt]=0.4;
 598   prodcutsval[2][ipt]=0.5;
 599   prodcutsval[3][ipt]=0.;
 600   prodcutsval[4][ipt]=0.;
 601   prodcutsval[5][ipt]=0.01;
 602   prodcutsval[6][ipt]=0.06;
 603   prodcutsval[7][ipt]=0.005;
 604   prodcutsval[8][ipt]=0.;
 605   prodcutsval[9][ipt]=0.;
 606   prodcutsval[10][ipt]=0.;
 607   prodcutsval[11][ipt]=0.05;
 608   prodcutsval[12][ipt]=0.4;
 609  }
 610  SetCuts(nvars,nptbins,prodcutsval);
 611
 612  AliAODPidHF* pidObjK=new AliAODPidHF();
 613  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
 614  pidObjK->SetSigma(sigmasK);
 615  pidObjK->SetAsym(kTRUE);
 616  pidObjK->SetMatch(1);
 617  pidObjK->SetTPC(kTRUE);
 618  pidObjK->SetTOF(kTRUE);
 619  pidObjK->SetITS(kTRUE);
 620  Double_t plimK[2]={0.5,0.8};
 621  pidObjK->SetPLimit(plimK,2);
 622  pidObjK->SetTOFdecide(kTRUE);
 623
 624  SetPidHF(pidObjK);
 625
 626  AliAODPidHF* pidObjpi=new AliAODPidHF();
 627  pidObjpi->SetTPC(kTRUE);
 628  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
 629  pidObjpi->SetSigma(sigmaspi);
 630  pidObjpi->SetTOFdecide(kTRUE);
 631  SetPidpion(pidObjpi);
 632
 633  AliAODPidHF* pidObjp=new AliAODPidHF();
 634  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
 635  pidObjp->SetSigma(sigmasp);
 636  pidObjp->SetAsym(kTRUE);
 637  pidObjp->SetMatch(1);
 638  pidObjp->SetTPC(kTRUE);
 639  pidObjp->SetTOF(kTRUE);
 640  pidObjp->SetITS(kTRUE);
 641  Double_t plimp[2]={1.,2.};
 642  pidObjp->SetPLimit(plimp,2);
 643  pidObjp->SetTOFdecide(kTRUE);
 644
 645  SetPidprot(pidObjp);
 646
 647  SetUsePID(kTRUE);
 648
 649  PrintAll();
 650
 651  for(Int_t iiv=0;iiv<nvars;iiv++){
 652   delete [] prodcutsval[iiv];
 653  }
 654  delete [] prodcutsval;
 655  prodcutsval=NULL;
 656  delete [] ptbins;
 657  ptbins=NULL;
 658
 659  delete pidObjp;
 660  pidObjp=NULL;
 661
 662  return;
 663 }
 664 //------------------
 665 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2010() {
 666
 667  SetName("LctopKpiProdCuts");
 668  SetTitle("Production cuts for Lc analysis");
 669
 670  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 671
 672  esdTrackCuts->SetRequireTPCRefit(kTRUE);
 673  esdTrackCuts->SetMinNClustersTPC(70);
 674  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 675                                           AliESDtrackCuts::kAny);
 676  esdTrackCuts->SetRequireITSRefit(kTRUE);
 677  esdTrackCuts->SetMinNClustersITS(4);
 678  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0100*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
 679  esdTrackCuts->SetEtaRange(-0.8,0.8);
 680  esdTrackCuts->SetMaxDCAToVertexXY(1.);
 681  esdTrackCuts->SetMaxDCAToVertexZ(1.);
 682  esdTrackCuts->SetPtRange(0.8,1.e10);
 683  AddTrackCuts(esdTrackCuts);
 684
 685  const Int_t nptbins=4;
 686  const Int_t nvars=13;
 687  Float_t* ptbins;
 688  ptbins=new Float_t[nptbins+1];
 689
 690  ptbins[0]=0.;
 691  ptbins[1]=2.;
 692  ptbins[2]=3.;
 693  ptbins[3]=4.;
 694  ptbins[4]=9999.;
 695
 696  SetGlobalIndex(nvars,nptbins);
 697  SetPtBins(nptbins+1,ptbins);
 698
 699  Float_t** prodcutsval;
 700  prodcutsval=new Float_t*[nvars];
 701  for(Int_t iv=0;iv<nvars;iv++){
 702   prodcutsval[iv]=new Float_t[nptbins];
 703  }
 704
 705  for(Int_t ipt=0;ipt<nptbins;ipt++){
 706   prodcutsval[0][ipt]=0.13;
 707   prodcutsval[1][ipt]=0.5;
 708   prodcutsval[2][ipt]=0.6;
 709   prodcutsval[3][ipt]=0.;
 710   prodcutsval[4][ipt]=0.;
 711   prodcutsval[5][ipt]=0.01;
 712   prodcutsval[6][ipt]=0.04;
 713   prodcutsval[7][ipt]=0.006;
 714   prodcutsval[8][ipt]=0.8;
 715   prodcutsval[9][ipt]=0.3;
 716   prodcutsval[10][ipt]=0.;
 717   prodcutsval[11][ipt]=0.05;
 718   prodcutsval[12][ipt]=0.4;
 719  }
 720  SetCuts(nvars,nptbins,prodcutsval);
 721
 722  AliAODPidHF* pidObjK=new AliAODPidHF();
 723  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
 724  pidObjK->SetSigma(sigmasK);
 725  pidObjK->SetAsym(kTRUE);
 726  pidObjK->SetMatch(1);
 727  pidObjK->SetTPC(kTRUE);
 728  pidObjK->SetTOF(kTRUE);
 729  pidObjK->SetITS(kTRUE);
 730  Double_t plimK[2]={0.5,0.8};
 731  pidObjK->SetPLimit(plimK,2);
 732
 733  SetPidHF(pidObjK);
 734
 735  AliAODPidHF* pidObjpi=new AliAODPidHF();
 736  pidObjpi->SetTPC(kTRUE);
 737  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
 738  pidObjpi->SetSigma(sigmaspi);
 739  SetPidpion(pidObjpi);
 740
 741  AliAODPidHF* pidObjp=new AliAODPidHF();
 742  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
 743  pidObjp->SetSigma(sigmasp);
 744  pidObjp->SetAsym(kTRUE);
 745  pidObjp->SetMatch(1);
 746  pidObjp->SetTPC(kTRUE);
 747  pidObjp->SetTOF(kTRUE);
 748  pidObjp->SetITS(kTRUE);
 749  Double_t plimp[2]={1.,2.};
 750  pidObjp->SetPLimit(plimp,2);
 751
 752  SetPidprot(pidObjp);
 753
 754  SetUsePID(kTRUE);
 755
 756  PrintAll();
 757
 758  for(Int_t iiv=0;iiv<nvars;iiv++){
 759   delete [] prodcutsval[iiv];
 760  }
 761  delete [] prodcutsval;
 762  prodcutsval=NULL;
 763  delete [] ptbins;
 764  ptbins=NULL;
 765
 766  delete pidObjp;
 767  pidObjp=NULL;
 768
 769  return;
 770 }
 771 //------------------
 772 AliKFParticle* AliRDHFCutsLctopKpi::ReconstructKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field,Bool_t constraint) const{
 773   // Method to construct the KF particle from the candidate
 774
 775  const Int_t nprongs=d->GetNProngs();
 776  if(nprongs<=0) return 0x0;
 777
 778  Int_t iprongs[nprongs];
 779  for(Int_t i=0;i<nprongs;i++) iprongs[i]=i;
 780
 781  Double_t mass[2]={0.,0.};
 782
 783  AliKFParticle *decay=d->ApplyVertexingKF(iprongs,nprongs,pdgs,constraint,field,mass);
 784  if(!decay) return 0x0;
 785  AliESDVertex *vertexESD = new AliESDVertex(decay->Parameters(),
 786                                             decay->CovarianceMatrix(),
 787                                             decay->GetChi2(),
 788                                             nprongs);
 789  Double_t pos[3],cov[6],chi2perNDF;
 790  vertexESD->GetXYZ(pos);
 791  vertexESD->GetCovMatrix(cov);
 792  chi2perNDF = vertexESD->GetChi2toNDF();
 793  delete vertexESD; vertexESD=NULL;
 794  AliAODVertex *vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 795  d->SetSecondaryVtx(vertexAOD);
 796  return decay;
 797 }
 798
 799 //------------------
 800 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2011() {
 801
 802   // Default 2010 PbPb cut object
 803   SetStandardCutsPbPb2010();
 804
 805   //
 806   // Enable all 2011 PbPb run triggers
 807   //
 808   SetTriggerClass("");
 809   ResetMaskAndEnableMBTrigger();
 810   EnableCentralTrigger();
 811   EnableSemiCentralTrigger();
 812 }
 813 //-----------------
 814
 815 Bool_t AliRDHFCutsLctopKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 816 {
 817   //
 818   //  // Checking if Dplus is in fiducial acceptance region
 819   //    //
 820   //
 821   if(pt > 5.) {
 822     // applying cut for pt > 5 GeV
 823    AliDebug(2,Form("pt of D+ = %f (> 5), cutting at |y| < 0.8",pt));
 824    if (TMath::Abs(y) > 0.8) return kFALSE;
 825
 826   } else {
 827    // appliying smooth cut for pt < 5 GeV
 828    Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 829    Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 830   AliDebug(2,Form("pt of D+ = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 831   if (y < minFiducialY || y > maxFiducialY) return kFALSE;
 832  }
 833   //
 834   return kTRUE;
 835 }
 836