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 #include <AliAnalysisManager.h>
  28 #include <AliInputEventHandler.h>
  29 #include <AliPIDResponse.h>
  30
  31 #include "AliRDHFCutsLctopKpi.h"
  32 #include "AliAODRecoDecayHF3Prong.h"
  33 #include "AliRDHFCuts.h"
  34 #include "AliAODTrack.h"
  35 #include "AliESDtrack.h"
  36 #include "AliKFParticle.h"
  37 #include "AliESDVertex.h"
  38
  39 using std::cout;
  40 using std::endl;
  41
  42 ClassImp(AliRDHFCutsLctopKpi)
  43
  44 //--------------------------------------------------------------------------
  45 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const char* name) :
  46 AliRDHFCuts(name),
  47 fPidObjprot(0),
  48 fPidObjpion(0),
  49 fUseImpParProdCorrCut(kFALSE),
  50 fPIDStrategy(kNSigma),
  51 fCutsStrategy(kStandard),
  52 fUseSpecialCut(kFALSE)
  53 {
  54   //
  55   // Default Constructor
  56   //
  57   Int_t nvars=13;
  58   SetNVars(nvars);
  59   TString varNames[13]={"inv. mass [GeV]",
  60                         "pTK [GeV/c]",
  61                         "pTP [GeV/c]",
  62                         "d0K [cm]   lower limit!",
  63                         "d0Pi [cm]  lower limit!",
  64                         "dist12 (cm)",
  65                         "sigmavert (cm)",
  66                         "dist prim-sec (cm)",
  67                         "pM=Max{pT1,pT2,pT3} (GeV/c)",
  68                         "cosThetaPoint",
  69                         "Sum d0^2 (cm^2)",
  70                         "dca cut (cm)",
  71                         "cut on pTpion [GeV/c]"};
  72   Bool_t isUpperCut[13]={kTRUE,
  73                          kFALSE,
  74                          kFALSE,
  75                          kFALSE,
  76                          kFALSE,
  77                          kFALSE,
  78                          kTRUE,
  79                          kFALSE,
  80                          kFALSE,
  81                          kFALSE,
  82                          kFALSE,
  83                          kTRUE,
  84                          kFALSE
  85                          };
  86   SetVarNames(nvars,varNames,isUpperCut);
  87   Bool_t forOpt[13]={kFALSE,
  88                      kTRUE,
  89                      kTRUE,
  90                      kFALSE,
  91                      kFALSE,
  92                      kFALSE,
  93                      kFALSE,
  94                      kTRUE,
  95                      kFALSE,
  96                      kFALSE,
  97                      kFALSE,
  98                      kFALSE,
  99                      kTRUE};
 100   SetVarsForOpt(4,forOpt);
 101   Float_t limits[2]={0,999999999.};
 102   SetPtBins(2,limits);
 103   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies)
 104       fPIDThreshold[ispecies]=0.;
 105 }
 106 //--------------------------------------------------------------------------
 107 AliRDHFCutsLctopKpi::AliRDHFCutsLctopKpi(const AliRDHFCutsLctopKpi &source) :
 108   AliRDHFCuts(source),
 109   fPidObjprot(0x0),
 110   fPidObjpion(0x0),
 111   fUseImpParProdCorrCut(source.fUseImpParProdCorrCut),
 112   fPIDStrategy(source.fPIDStrategy),
 113   fCutsStrategy(source.fCutsStrategy),
 114   fUseSpecialCut(source.fUseSpecialCut)
 115 {
 116   //
 117   // Copy constructor
 118   //
 119   if (source.fPidObjprot) fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
 120   else fPidObjprot = new AliAODPidHF();
 121   if (source.fPidObjpion) fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
 122   else fPidObjpion = new AliAODPidHF();
 123   memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
 124 }
 125 //--------------------------------------------------------------------------
 126 AliRDHFCutsLctopKpi &AliRDHFCutsLctopKpi::operator=(const AliRDHFCutsLctopKpi &source)
 127 {
 128   //
 129   // assignment operator
 130   //
 131   if(this != &source) {
 132
 133     AliRDHFCuts::operator=(source);
 134     delete fPidObjprot;
 135     fPidObjprot = new AliAODPidHF(*(source.fPidObjprot));
 136     delete fPidObjpion;
 137     fPidObjpion = new AliAODPidHF(*(source.fPidObjpion));
 138     fPIDStrategy=source.fPIDStrategy;
 139     fCutsStrategy=source.fCutsStrategy;
 140     memcpy(fPIDThreshold,source.fPIDThreshold,AliPID::kSPECIES*sizeof(Double_t));
 141   }
 142
 143   return *this;
 144 }
 145 //---------------------------------------------------------------------------
 146 AliRDHFCutsLctopKpi::~AliRDHFCutsLctopKpi() {
 147  //
 148  //  // Default Destructor
 149  //
 150  if(fPidObjpion){
 151   delete fPidObjpion;
 152   fPidObjpion=0;
 153  }
 154  if(fPidObjprot){
 155   delete fPidObjprot;
 156   fPidObjprot=0;
 157  }
 158
 159 }
 160
 161 //---------------------------------------------------------------------------
 162 void AliRDHFCutsLctopKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters, AliAODEvent *aod) {
 163   //
 164   // Fills in vars the values of the variables
 165   //
 166
 167   if(nvars!=fnVarsForOpt) {
 168     printf("AliRDHFCutsLctopKpi::GetCutsVarsForOpt: wrong number of variables\n");
 169     return;
 170   }
 171
 172   AliAODRecoDecayHF3Prong *dd = (AliAODRecoDecayHF3Prong*)d;
 173
 174   Int_t iter=-1;
 175   if(fVarsForOpt[0]){
 176     iter++;
 177     vars[iter]=dd->InvMassLcpKpi();
 178   }
 179   if(fVarsForOpt[1]){
 180     iter++;
 181     for(Int_t iprong=0;iprong<3;iprong++){
 182       if(TMath::Abs(pdgdaughters[iprong])==321) {
 183         vars[iter]=dd->PtProng(iprong);
 184       }
 185     }
 186   }
 187   if(fVarsForOpt[2]){
 188     iter++;
 189     for(Int_t iprong=0;iprong<3;iprong++){
 190       if(TMath::Abs(pdgdaughters[iprong])==2212) {
 191         vars[iter]=dd->PtProng(iprong);
 192       }
 193     }
 194   }
 195   if(fVarsForOpt[3]){
 196     iter++;
 197     for(Int_t iprong=0;iprong<3;iprong++){
 198       if(TMath::Abs(pdgdaughters[iprong])==2212) {
 199         vars[iter]=dd->Getd0Prong(iprong);
 200       }
 201     }
 202   }
 203   if(fVarsForOpt[4]){
 204     iter++;
 205     for(Int_t iprong=0;iprong<3;iprong++){
 206       if(TMath::Abs(pdgdaughters[iprong])==211) {
 207         vars[iter]=dd->Getd0Prong(iprong);
 208       }
 209     }
 210   }
 211   if(fVarsForOpt[5]){
 212     iter++;
 213     vars[iter]=dd->GetDist12toPrim();
 214   }
 215   if(fVarsForOpt[6]){
 216     iter++;
 217     vars[iter]=dd->GetSigmaVert(aod);
 218   }
 219   if(fVarsForOpt[7]){
 220     iter++;
 221     vars[iter] = dd->DecayLength();
 222   }
 223   if(fVarsForOpt[8]){
 224     iter++;
 225     Float_t ptmax=0;
 226     for(Int_t i=0;i<3;i++){
 227       if(dd->PtProng(i)>ptmax)ptmax=dd->PtProng(i);
 228     }
 229     vars[iter]=ptmax;
 230   }
 231   if(fVarsForOpt[9]){
 232     iter++;
 233     vars[iter]=dd->CosPointingAngle();
 234   }
 235   if(fVarsForOpt[10]){
 236     iter++;
 237     vars[iter]=dd->Getd0Prong(0)*dd->Getd0Prong(0)+dd->Getd0Prong(1)*dd->Getd0Prong(1)+dd->Getd0Prong(2)*dd->Getd0Prong(2);
 238   }
 239   if(fVarsForOpt[11]){
 240     iter++;
 241     vars[iter]=dd->GetDCA();
 242   }
 243   if(fVarsForOpt[12]){
 244     iter++;
 245     for(Int_t iprong=0;iprong<3;iprong++){
 246       if(TMath::Abs(pdgdaughters[iprong])==211) {
 247         vars[iter]=dd->PtProng(iprong);
 248       }
 249     }
 250   }
 251
 252   return;
 253 }
 254 //---------------------------------------------------------------------------
 255 Int_t AliRDHFCutsLctopKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEvent *aod) {
 256   //
 257   // Apply selection
 258   //
 259
 260   if(!fCutsRD){
 261     AliError("Cut matrice not inizialized. Exit...\n");
 262     return 0;
 263   }
 264   //PrintAll();
 265   AliAODRecoDecayHF3Prong* d=(AliAODRecoDecayHF3Prong*)obj;
 266
 267   if(!d){
 268     AliError("AliAODRecoDecayHF3Prong null \n");
 269     return 0;
 270   }
 271
 272
 273   if(fKeepSignalMC) if(IsSignalMC(d,aod,4122)) return 3;
 274
 275   Int_t returnvalue=3;
 276   Int_t returnvaluePID=3;
 277
 278   if(d->Pt()<fMinPtCand) return 0;
 279   if(d->Pt()>fMaxPtCand) return 0;
 280
 281   if(fUseTrackSelectionWithFilterBits && d->HasBadDaughters()) return 0;
 282
 283
 284   // selection on daughter tracks
 285   if(selectionLevel==AliRDHFCuts::kAll ||
 286      selectionLevel==AliRDHFCuts::kTracks) {
 287     if(!AreDaughtersSelected(d)) return 0;
 288   }
 289
 290
 291   // PID selection
 292   if(selectionLevel==AliRDHFCuts::kAll ||
 293      selectionLevel==AliRDHFCuts::kCandidate||
 294      selectionLevel==AliRDHFCuts::kPID) {
 295     switch (fPIDStrategy) {
 296     case kNSigma:
 297       returnvaluePID = IsSelectedPID(d);
 298       break;
 299      case kNSigmaPbPb:
 300       returnvaluePID = IsSelectedNSigmaPbPb(d);
 301       break;
 302     case kCombined:
 303       returnvaluePID = IsSelectedCombinedPID(d);
 304       break;
 305     case kCombinedSoft:
 306       returnvaluePID = IsSelectedCombinedPIDSoft(d);
 307       break;
 308     case kNSigmaStrong:
 309       returnvaluePID = IsSelectedPIDStrong(d);
 310       break;
 311     case kCombinedpPb:
 312       returnvaluePID = IsSelectedCombinedPIDpPb(d);
 313       break;
 314     case kCombinedpPb2:
 315       returnvaluePID = IsSelectedCombinedPIDpPb2(d);
 316       break;
 317     }
 318     fIsSelectedPID=returnvaluePID;
 319   }
 320   //  if(fUsePID || selectionLevel==AliRDHFCuts::kPID) returnvaluePID = IsSelectedCombinedPID(d);   // to test!!
 321   if(returnvaluePID==0) return 0;
 322
 323
 324
 325
 326   // selection on candidate
 327   if(selectionLevel==AliRDHFCuts::kAll ||
 328      selectionLevel==AliRDHFCuts::kCandidate) {
 329
 330     Double_t pt=d->Pt();
 331
 332     Int_t ptbin=PtBin(pt);
 333
 334     Double_t mLcpKpi=0.,mLcpiKp=0.;
 335     Int_t okLcpKpi=1,okLcpiKp=1;
 336
 337     Double_t mLcPDG = TDatabasePDG::Instance()->GetParticle(4122)->Mass();
 338
 339     mLcpKpi=d->InvMassLcpKpi();
 340     mLcpiKp=d->InvMassLcpiKp();
 341
 342     if(TMath::Abs(mLcpKpi-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpKpi = 0;
 343     if(TMath::Abs(mLcpiKp-mLcPDG)>fCutsRD[GetGlobalIndex(0,ptbin)]) okLcpiKp = 0;
 344     if(!okLcpKpi && !okLcpiKp) return 0;
 345
 346   switch (fCutsStrategy) {
 347
 348     case kStandard:
 349     if(TMath::Abs(d->PtProng(1)) < fCutsRD[GetGlobalIndex(1,ptbin)] || TMath::Abs(d->Getd0Prong(1))<fCutsRD[GetGlobalIndex(3,ptbin)]) return 0;//Kaon
 350     if(d->Pt()>=3. && d->PProng(1)<0.55) return 0;
 351     if(fUseSpecialCut) {
 352       if(TMath::Abs(d->PtProng(0)) < TMath::Abs(d->PtProng(2)) )okLcpKpi=0;
 353       if(TMath::Abs(d->PtProng(2)) < TMath::Abs(d->PtProng(0)) )okLcpiKp=0;
 354     }
 355     if((TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(12,ptbin)])) okLcpKpi=0;
 356     if((TMath::Abs(d->PtProng(2)) < fCutsRD[GetGlobalIndex(2,ptbin)]) || (TMath::Abs(d->PtProng(0)) < fCutsRD[GetGlobalIndex(12,ptbin)]))okLcpiKp=0;
 357     if(!okLcpKpi && !okLcpiKp) return 0;
 358     //2track cuts
 359     if(d->GetDist12toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]|| d->GetDist23toPrim()<fCutsRD[GetGlobalIndex(5,ptbin)]) return 0;
 360     if(d->GetDist12toPrim()>0.5) return 0;
 361     if(d->GetDist23toPrim()>0.5) return 0;
 362     if(fUseImpParProdCorrCut){
 363       if(d->Getd0Prong(0)*d->Getd0Prong(1)<0. && d->Getd0Prong(2)*d->Getd0Prong(1)<0.) return 0;
 364     }
 365     //sec vert
 366     if(d->DecayLength()<fCutsRD[GetGlobalIndex(7,ptbin)]) return 0;
 367     if(d->DecayLength()>0.5) return 0;
 368
 369   //  Double_t sumd0s=d->Getd0Prong(0)*d->Getd0Prong(0)+d->Getd0Prong(1)*d->Getd0Prong(1)+d->Getd0Prong(2)*d->Getd0Prong(2);
 370   //  if(sumd0s<fCutsRD[GetGlobalIndex(10,ptbin)]) return 0;
 371     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;
 372
 373     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;
 374     if(d->CosPointingAngle()< fCutsRD[GetGlobalIndex(9,ptbin)]) return 0;
 375     if(d->GetSigmaVert(aod)>fCutsRD[GetGlobalIndex(6,ptbin)]) return 0;
 376
 377     //DCA
 378     for(Int_t i=0;i<3;i++) if(d->GetDCA(i)>fCutsRD[GetGlobalIndex(11,ptbin)]) return 0;
 379
 380     break;
 381
 382    case kKF:
 383     Int_t pdgs[3]={0,321,0};
 384     Bool_t constraint=kFALSE;
 385     if(fCutsRD[GetGlobalIndex(1,ptbin)]>0.) constraint=kTRUE;
 386     Double_t field=aod->GetMagneticField();
 387     if (returnvaluePID==1 || returnvaluePID==3){
 388
 389       pdgs[0]=2122;pdgs[2]=211;
 390       AliKFParticle *lc1=ReconstructKF(d,pdgs,field,constraint);
 391       if(!lc1){
 392         okLcpKpi=0;
 393       }else{
 394         if(lc1->GetChi2()/lc1->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)]) okLcpKpi=0;
 395       }
 396     } else if(returnvaluePID>=2){
 397
 398       pdgs[0]=211;pdgs[2]=2212;
 399       AliKFParticle *lc2=ReconstructKF(d,pdgs,field,constraint);
 400       if(!lc2){
 401         okLcpiKp=0;
 402       }else{
 403         if(lc2->GetChi2()/lc2->GetNDF()>fCutsRD[GetGlobalIndex(2,ptbin)])okLcpiKp=0;
 404       }
 405     }
 406     break;
 407
 408    }
 409
 410     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 411     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 412     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 413
 414   }
 415
 416
 417   Int_t returnvalueTot=CombinePIDCuts(returnvalue,returnvaluePID);
 418   return returnvalueTot;
 419 }
 420 //---------------------------------------------------------------------------
 421 Int_t AliRDHFCutsLctopKpi::IsSelectedPID(AliAODRecoDecayHF* obj) {
 422
 423
 424     if(!fUsePID || !obj) return 3;
 425     Int_t okLcpKpi=0,okLcpiKp=0;
 426     Int_t returnvalue=0;
 427     Bool_t isPeriodd=fPidHF->GetOnePad();
 428     Bool_t isMC=fPidHF->GetMC();
 429     Bool_t ispion0=kTRUE,ispion2=kTRUE;
 430     Bool_t isproton0=kFALSE,isproton2=kFALSE;
 431     Bool_t iskaon1=kFALSE;
 432     if(isPeriodd) {
 433      fPidObjprot->SetOnePad(kTRUE);
 434      fPidObjpion->SetOnePad(kTRUE);
 435     }
 436     if(isMC) {
 437      fPidObjprot->SetMC(kTRUE);
 438      fPidObjpion->SetMC(kTRUE);
 439     }
 440
 441    if(fPidObjprot->GetPidResponse()==0x0){
 442       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 443       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 444       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 445       fPidObjprot->SetPidResponse(pidResp);
 446     }
 447     if(fPidObjpion->GetPidResponse()==0x0){
 448       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 449       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 450       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 451       fPidObjpion->SetPidResponse(pidResp);
 452     }
 453     if(fPidHF->GetPidResponse()==0x0){
 454       AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 455       AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 456       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 457       fPidHF->SetPidResponse(pidResp);
 458     }
 459
 460     for(Int_t i=0;i<3;i++){
 461      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
 462      if(!track) return 0;
 463      // identify kaon
 464      if(track->P()<0.55){
 465       fPidHF->SetTOF(kFALSE);
 466       fPidHF->SetTOFdecide(kFALSE);
 467      }
 468      if(i==1) {
 469       Int_t isKaon=fPidHF->MakeRawPid(track,3);
 470       if(isKaon>=1) iskaon1=kTRUE;
 471       if(track->P()<0.55){
 472        fPidHF->SetTOF(kTRUE);
 473        fPidHF->SetTOFdecide(kTRUE);
 474       }
 475
 476       if(!iskaon1) return 0;
 477
 478      }else{
 479      //pion or proton
 480     if(track->P()<1.){
 481       fPidObjprot->SetTOF(kFALSE);
 482       fPidObjprot->SetTOFdecide(kFALSE);
 483      }
 484
 485      Int_t isProton=fPidObjprot->MakeRawPid(track,4);
 486
 487
 488      Int_t isPion=fPidObjpion->MakeRawPid(track,2);
 489
 490      if(track->P()<1.){
 491       fPidObjprot->SetTOF(kTRUE);
 492       fPidObjprot->SetTOFdecide(kTRUE);
 493      }
 494
 495
 496      if(i==0) {
 497       if(isPion<0) ispion0=kFALSE;
 498       if(isProton>=1) isproton0=kTRUE;
 499
 500      }
 501       if(!ispion0 && !isproton0) return 0;
 502      if(i==2) {
 503       if(isPion<0) ispion2=kFALSE;
 504       if(isProton>=1) isproton2=kTRUE;
 505      }
 506
 507     }
 508    }
 509
 510     if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
 511     if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
 512     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 513     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 514     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 515
 516  return returnvalue;
 517 }
 518
 519 //--------------------------------------------------------------------------
 520
 521 Int_t AliRDHFCutsLctopKpi::IsSelectedNSigmaPbPb(AliAODRecoDecayHF* obj) {
 522
 523
 524   if(!fUsePID || !obj) return 3;
 525   Int_t okLcpKpi=0,okLcpiKp=0;
 526   Int_t returnvalue=0;
 527   Bool_t isPeriodd=fPidHF->GetOnePad();
 528   Bool_t isMC=fPidHF->GetMC();
 529   Bool_t ispion0=kTRUE,ispion2=kTRUE;
 530   Bool_t isproton0=kFALSE,isproton2=kFALSE;
 531   Bool_t iskaon1=kFALSE;
 532   if(isPeriodd) {
 533     fPidObjprot->SetOnePad(kTRUE);
 534     fPidObjpion->SetOnePad(kTRUE);
 535   }
 536   if(isMC) {
 537     fPidObjprot->SetMC(kTRUE);
 538     fPidObjpion->SetMC(kTRUE);
 539   }
 540
 541   if(fPidObjprot->GetPidResponse()==0x0){
 542     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 543     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 544     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 545     fPidObjprot->SetPidResponse(pidResp);
 546   }
 547   if(fPidObjpion->GetPidResponse()==0x0){
 548     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 549     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 550     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 551     fPidObjpion->SetPidResponse(pidResp);
 552   }
 553   if(fPidHF->GetPidResponse()==0x0){
 554     AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
 555     AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
 556     AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
 557     fPidHF->SetPidResponse(pidResp);
 558   }
 559
 560   for(Int_t i=0;i<3;i++){
 561     AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
 562     if(!track) return 0;
 563
 564     if(i==1) {
 565       //kaon
 566       Int_t isKaon=fPidHF->MakeRawPid(track,3);
 567       if(isKaon>=1) iskaon1=kTRUE;
 568       if(!iskaon1) return 0;
 569     }
 570
 571     else {
 572       //pion or proton
 573       Int_t isProton=fPidObjprot->MakeRawPid(track,4);
 574       Int_t isPion=fPidObjpion->MakeRawPid(track,2);
 575
 576       if(i==0) {
 577         if(isPion<0) ispion0=kFALSE;
 578         if(isProton>=1) isproton0=kTRUE;
 579       }
 580       if(!ispion0 && !isproton0) return 0;
 581
 582       if(i==2) {
 583         if(isPion<0) ispion2=kFALSE;
 584         if(isProton>=1) isproton2=kTRUE;
 585       }
 586     }
 587   }
 588
 589   if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
 590   if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
 591   if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 592   if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 593   if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 594
 595   return returnvalue;
 596 }
 597
 598
 599
 600
 601 //---------------------------------------------------------------------------
 602 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPID(AliAODRecoDecayHF* obj) {
 603
 604     if(!fUsePID || !obj) {return 3;}
 605     Int_t okLcpKpi=0,okLcpiKp=0;
 606     Int_t returnvalue=0;
 607     Bool_t isPeriodd=fPidHF->GetOnePad();
 608     Bool_t isMC=fPidHF->GetMC();
 609
 610     if(isPeriodd) {
 611             fPidObjprot->SetOnePad(kTRUE);
 612             fPidObjpion->SetOnePad(kTRUE);
 613     }
 614     if(isMC) {
 615             fPidObjprot->SetMC(kTRUE);
 616             fPidObjpion->SetMC(kTRUE);
 617     }
 618
 619     AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
 620     AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
 621     AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
 622     if (!track0 || !track1 || !track2) return 0;
 623     Double_t prob0[AliPID::kSPECIES];
 624     Double_t prob1[AliPID::kSPECIES];
 625     Double_t prob2[AliPID::kSPECIES];
 626     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 627     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
 628     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 629
 630    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 631     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
 632    if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 633
 634     if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 635     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
 636    if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 637
 638     if(fPIDThreshold[AliPID::kPion]>0. && fPIDThreshold[AliPID::kKaon]>0. && fPIDThreshold[AliPID::kProton]>0.){
 639     okLcpiKp=  (prob0[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ])
 640              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
 641              &&(prob2[AliPID::kProton]>fPIDThreshold[AliPID::kProton]);
 642     okLcpKpi=  (prob0[AliPID::kProton]>fPIDThreshold[AliPID::kProton])
 643              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
 644              &&(prob2[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ]);
 645    }else{
 646                     //pion or proton
 647
 648
 649     if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
 650     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;
 651     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1;
 652             }
 653            }
 654
 655     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
 656     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
 657     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
 658
 659     return returnvalue;
 660 }
 661 //-----------------------
 662 Int_t AliRDHFCutsLctopKpi::CombinePIDCuts(Int_t returnvalue, Int_t returnvaluePID) const {
 663
 664  Int_t returnvalueTot=0;
 665  Int_t okLcpKpi=0,okLcpiKp=0;
 666  if(returnvaluePID==1){
 667    if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
 668  }
 669  if(returnvaluePID==2){
 670    if(returnvalue>=2) okLcpiKp=1;
 671  }
 672  if(returnvaluePID==3 && returnvalue>0){
 673   if(returnvalue==1 || returnvalue==3) okLcpKpi=1;
 674   if(returnvalue>=2) okLcpiKp=1;
 675  }
 676
 677  if(okLcpKpi) returnvalueTot=1; //cuts passed as Lc->pKpi
 678  if(okLcpiKp) returnvalueTot=2; //cuts passed as Lc->piKp
 679  if(okLcpKpi && okLcpiKp) returnvalueTot=3; //cuts passed as both pKpi and piKp
 680  return returnvalueTot;
 681 }
 682 //----------------------------------
 683 void AliRDHFCutsLctopKpi::SetStandardCutsPP2010() {
 684
 685  SetName("LctopKpiProdCuts");
 686  SetTitle("Production cuts for Lc analysis");
 687
 688  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 689  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
 690  esdTrackCuts->SetRequireTPCRefit(kTRUE);
 691  esdTrackCuts->SetMinNClustersTPC(70);
 692  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 693                                           AliESDtrackCuts::kAny);
 694  esdTrackCuts->SetRequireITSRefit(kTRUE);
 695  esdTrackCuts->SetMinNClustersITS(4);
 696  esdTrackCuts->SetMinDCAToVertexXY(0.);
 697  esdTrackCuts->SetEtaRange(-0.8,0.8);
 698  esdTrackCuts->SetPtRange(0.3,1.e10);
 699  AddTrackCuts(esdTrackCuts);
 700  delete esdTrackCuts;
 701  esdTrackCuts=NULL;
 702
 703  const Int_t nptbins=8;
 704  const Int_t nvars=13;
 705  Float_t* ptbins;
 706  ptbins=new Float_t[nptbins+1];
 707
 708  ptbins[0]=1.;
 709  ptbins[1]=2.;
 710  ptbins[2]=3.;
 711  ptbins[3]=4.;
 712  ptbins[4]=5.;
 713  ptbins[5]=6.;
 714  ptbins[6]=8.;
 715  ptbins[7]=10.;
 716  ptbins[8]=20.;
 717
 718  SetGlobalIndex(nvars,nptbins);
 719  SetPtBins(nptbins+1,ptbins);
 720
 721  Float_t** prodcutsval;
 722  prodcutsval=new Float_t*[nvars];
 723  for(Int_t iv=0;iv<nvars;iv++){
 724   prodcutsval[iv]=new Float_t[nptbins];
 725  }
 726
 727  for(Int_t ipt=0;ipt<nptbins;ipt++){
 728   prodcutsval[0][ipt]=0.13;
 729   prodcutsval[1][ipt]=0.4;
 730   prodcutsval[2][ipt]=0.4;
 731   prodcutsval[3][ipt]=0.;
 732   prodcutsval[4][ipt]=0.;
 733   prodcutsval[5][ipt]=0.;
 734   prodcutsval[6][ipt]=0.06;
 735   prodcutsval[7][ipt]=0.005;
 736   prodcutsval[8][ipt]=0.;
 737   prodcutsval[9][ipt]=0.;
 738   prodcutsval[10][ipt]=0.;
 739   prodcutsval[11][ipt]=0.05;
 740   prodcutsval[12][ipt]=0.4;
 741  }
 742  SetCuts(nvars,nptbins,prodcutsval);
 743
 744  AliAODPidHF* pidObjK=new AliAODPidHF();
 745  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
 746  pidObjK->SetSigma(sigmasK);
 747  pidObjK->SetAsym(kTRUE);
 748  pidObjK->SetMatch(1);
 749  pidObjK->SetTPC(kTRUE);
 750  pidObjK->SetTOF(kTRUE);
 751  pidObjK->SetITS(kTRUE);
 752  Double_t plimK[2]={0.5,0.8};
 753  pidObjK->SetPLimit(plimK,2);
 754  pidObjK->SetTOFdecide(kTRUE);
 755
 756  SetPidHF(pidObjK);
 757
 758  AliAODPidHF* pidObjpi=new AliAODPidHF();
 759  pidObjpi->SetTPC(kTRUE);
 760  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
 761  pidObjpi->SetSigma(sigmaspi);
 762  pidObjpi->SetTOFdecide(kTRUE);
 763  SetPidpion(pidObjpi);
 764
 765  AliAODPidHF* pidObjp=new AliAODPidHF();
 766  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
 767  pidObjp->SetSigma(sigmasp);
 768  pidObjp->SetAsym(kTRUE);
 769  pidObjp->SetMatch(1);
 770  pidObjp->SetTPC(kTRUE);
 771  pidObjp->SetTOF(kTRUE);
 772  pidObjp->SetITS(kTRUE);
 773  Double_t plimp[2]={1.,2.};
 774  pidObjp->SetPLimit(plimp,2);
 775  pidObjp->SetTOFdecide(kTRUE);
 776
 777  SetPidprot(pidObjp);
 778
 779  SetUsePID(kTRUE);
 780  SetOptPileup(kTRUE);
 781
 782  PrintAll();
 783
 784  for(Int_t iiv=0;iiv<nvars;iiv++){
 785   delete [] prodcutsval[iiv];
 786  }
 787  delete [] prodcutsval;
 788  prodcutsval=NULL;
 789  delete [] ptbins;
 790  ptbins=NULL;
 791
 792  delete pidObjK;
 793  pidObjK=NULL;
 794  delete pidObjpi;
 795  pidObjpi=NULL;
 796  delete pidObjp;
 797  pidObjp=NULL;
 798
 799  return;
 800 }
 801 //------------------
 802 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2010() {
 803
 804  SetName("LctopKpiProdCuts");
 805  SetTitle("Production cuts for Lc analysis");
 806
 807  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
 808
 809  esdTrackCuts->SetRequireTPCRefit(kTRUE);
 810  esdTrackCuts->SetMinNClustersTPC(70);
 811  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
 812                                           AliESDtrackCuts::kAny);
 813  esdTrackCuts->SetRequireITSRefit(kTRUE);
 814  esdTrackCuts->SetMinNClustersITS(4);
 815  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0100*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");
 816  esdTrackCuts->SetEtaRange(-0.8,0.8);
 817  esdTrackCuts->SetMaxDCAToVertexXY(1.);
 818  esdTrackCuts->SetMaxDCAToVertexZ(1.);
 819  esdTrackCuts->SetPtRange(0.49,1.e10);
 820  AddTrackCuts(esdTrackCuts);
 821  delete esdTrackCuts;
 822  esdTrackCuts=NULL;
 823
 824  const Int_t nptbins=8;
 825  const Int_t nvars=13;
 826  Float_t* ptbins;
 827  ptbins=new Float_t[nptbins+1];
 828
 829  ptbins[0]=1.;
 830  ptbins[1]=2.;
 831  ptbins[2]=3.;
 832  ptbins[3]=4.;
 833  ptbins[4]=5.;
 834  ptbins[5]=6.;
 835  ptbins[6]=8.;
 836  ptbins[7]=10.;
 837  ptbins[8]=99999.;
 838
 839
 840  SetGlobalIndex(nvars,nptbins);
 841  SetPtBins(nptbins+1,ptbins);
 842
 843  Float_t** prodcutsval;
 844  prodcutsval=new Float_t*[nvars];
 845  for(Int_t iv=0;iv<nvars;iv++){
 846   prodcutsval[iv]=new Float_t[nptbins];
 847  }
 848
 849  for(Int_t ipt=0;ipt<nptbins;ipt++){
 850   prodcutsval[0][ipt]=0.13;
 851   prodcutsval[1][ipt]=0.5;
 852   prodcutsval[2][ipt]=0.6;
 853   prodcutsval[3][ipt]=0.;
 854   prodcutsval[4][ipt]=0.;
 855   prodcutsval[5][ipt]=0.01;
 856   prodcutsval[6][ipt]=0.04;
 857   prodcutsval[7][ipt]=0.006;
 858   prodcutsval[8][ipt]=0.8;
 859   prodcutsval[9][ipt]=0.3;
 860   prodcutsval[10][ipt]=0.;
 861   prodcutsval[11][ipt]=0.05;
 862   prodcutsval[12][ipt]=0.4;
 863  }
 864  SetCuts(nvars,nptbins,prodcutsval);
 865
 866  AliAODPidHF* pidObj=new AliAODPidHF();
 867  pidObj->SetTPC(kTRUE);
 868  pidObj->SetTOF(kTRUE);
 869  SetPidHF(pidObj);
 870  SetPidpion(pidObj);
 871  SetPidprot(pidObj);
 872
 873
 874  // bayesian pid
 875  GetPidHF()->SetUseCombined(kTRUE);
 876  GetPidHF()->SetUseDefaultPriors(kTRUE);
 877  GetPidHF()->SetCombDetectors(AliAODPidHF::kTPCTOF);
 878  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies){
 879    SetPIDThreshold(static_cast<AliPID::EParticleType>(ispecies),0);
 880  }
 881  SetPIDStrategy(AliRDHFCutsLctopKpi::kCombinedpPb);
 882  SetUsePID(kTRUE);
 883
 884
 885  PrintAll();
 886
 887  for(Int_t iiv=0;iiv<nvars;iiv++){
 888   delete [] prodcutsval[iiv];
 889  }
 890  delete [] prodcutsval;
 891  prodcutsval=NULL;
 892  delete [] ptbins;
 893  ptbins=NULL;
 894
 895  delete pidObj;
 896  pidObj=NULL;
 897
 898  return;
 899 }
 900 //------------------
 901 AliKFParticle* AliRDHFCutsLctopKpi::ReconstructKF(AliAODRecoDecayHF3Prong *d,Int_t *pdgs,Double_t field,Bool_t constraint) const{
 902   // Method to construct the KF particle from the candidate
 903
 904  const Int_t nprongs=d->GetNProngs();
 905  if(nprongs<=0) return 0x0;
 906
 907  Int_t iprongs[nprongs];
 908  for(Int_t i=0;i<nprongs;i++) iprongs[i]=i;
 909
 910  Double_t mass[2]={0.,0.};
 911
 912  AliKFParticle *decay=d->ApplyVertexingKF(iprongs,nprongs,pdgs,constraint,field,mass);
 913  if(!decay) return 0x0;
 914  AliESDVertex *vertexESD = new AliESDVertex(decay->Parameters(),
 915                                             decay->CovarianceMatrix(),
 916                                             decay->GetChi2(),
 917                                             nprongs);
 918  Double_t pos[3],cov[6],chi2perNDF;
 919  vertexESD->GetXYZ(pos);
 920  vertexESD->GetCovMatrix(cov);
 921  chi2perNDF = vertexESD->GetChi2toNDF();
 922  delete vertexESD; vertexESD=NULL;
 923  AliAODVertex *vertexAOD = new AliAODVertex(pos,cov,chi2perNDF,0x0,-1,AliAODVertex::kUndef,nprongs);
 924  d->SetSecondaryVtx(vertexAOD);
 925  return decay;
 926 }
 927
 928 //------------------
 929 void AliRDHFCutsLctopKpi::SetStandardCutsPbPb2011() {
 930
 931   // Default 2010 PbPb cut object
 932   SetStandardCutsPbPb2010();
 933
 934   //
 935   // Enable all 2011 PbPb run triggers
 936   //
 937   SetTriggerClass("");
 938   ResetMaskAndEnableMBTrigger();
 939   EnableCentralTrigger();
 940   EnableSemiCentralTrigger();
 941 }
 942 //-----------------
 943
 944 Bool_t AliRDHFCutsLctopKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 945 {
 946   //
 947   //  // Checking if Lc is in fiducial acceptance region
 948   //    //
 949   //
 950   if(fMaxRapidityCand>-998.){
 951     if(TMath::Abs(y) > fMaxRapidityCand) return kFALSE;
 952     else return kTRUE;
 953   }
 954
 955   if(pt > 5.) {
 956     // applying cut for pt > 5 GeV
 957    AliDebug(2,Form("pt of Lc = %f (> 5), cutting at |y| < 0.8",pt));
 958    if (TMath::Abs(y) > 0.8) return kFALSE;
 959
 960   } else {
 961    // appliying smooth cut for pt < 5 GeV
 962    Double_t maxFiducialY = -0.2/15*pt*pt+1.9/15*pt+0.5;
 963    Double_t minFiducialY = 0.2/15*pt*pt-1.9/15*pt-0.5;
 964   AliDebug(2,Form("pt of Lc = %f (< 5), cutting  according to the fiducial zone [%f, %f]\n",pt,minFiducialY,maxFiducialY));
 965   if (y < minFiducialY || y > maxFiducialY) return kFALSE;
 966  }
 967   //
 968   return kTRUE;
 969 }
 970 //--------------------------------------------------------
 971 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDSoft(AliAODRecoDecayHF* obj) {
 972  if(!fUsePID || !obj) {return 3;}
 973  Int_t okLcpKpi=0,okLcpiKp=0;
 974  Int_t returnvalue=0;
 975
 976  AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
 977  AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
 978  AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
 979  if (!track0 || !track1 || !track2) return 0;
 980  Double_t prob0[AliPID::kSPECIES];
 981  Double_t prob1[AliPID::kSPECIES];
 982  Double_t prob2[AliPID::kSPECIES];
 983
 984  Bool_t isTOF0=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(0));
 985  Bool_t isTOF1=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(1));
 986  Bool_t isTOF2=fPidHF->CheckTOFPIDStatus((AliAODTrack*)obj->GetDaughter(2));
 987
 988 Bool_t isK1=kFALSE;
 989  if(isTOF1){ //kaon
 990   if(track1->P()<1.8) {
 991     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
 992     if(obj->Pt()<3. && track1->P()<0.55) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
 993     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
 994
 995   }else{
 996     AliAODTrack *trackaod1=(AliAODTrack*)(obj->GetDaughter(1));
 997     if(trackaod1->P()<0.55){
 998       fPidHF->SetTOF(kFALSE);
 999       fPidHF->SetTOFdecide(kFALSE);
1000      }
1001     Int_t isKaon=fPidHF->MakeRawPid(trackaod1,3);
1002     if(isKaon>=1) isK1=kTRUE;
1003     if(trackaod1->P()<0.55){
1004       fPidHF->SetTOF(kTRUE);
1005       fPidHF->SetTOFdecide(kTRUE);
1006      }
1007   }
1008  }else{
1009   if(track1->P()<0.8){
1010     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1011     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob0);
1012   }else{
1013     AliAODTrack *trackaod1=(AliAODTrack*)(obj->GetDaughter(1));
1014      if(trackaod1->P()<0.55){
1015       fPidHF->SetTOF(kFALSE);
1016       fPidHF->SetTOFdecide(kFALSE);
1017      }
1018     Int_t isKaon=fPidHF->MakeRawPid(trackaod1,3);
1019     if(isKaon>=1) isK1=kTRUE;
1020      if(trackaod1->P()<0.55){
1021       fPidHF->SetTOF(kTRUE);
1022       fPidHF->SetTOFdecide(kTRUE);
1023      }
1024   }
1025  }
1026
1027  Bool_t ispi0=kFALSE;
1028  Bool_t isp0=kFALSE;
1029  Bool_t ispi2=kFALSE;
1030  Bool_t isp2=kFALSE;
1031
1032  if(isTOF0){ //proton
1033   if(track0->P()<2.2) {
1034     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1035     if(obj->Pt()<3. && track0->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1036     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1037   }else{
1038    AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(0));
1039    if(trackaod0->P()<1.){
1040       fPidObjprot->SetTOF(kFALSE);
1041       fPidObjprot->SetTOFdecide(kFALSE);
1042    }
1043    Int_t isProton=fPidObjprot->MakeRawPid(trackaod0,4);
1044    if(isProton>=1) isp0=kTRUE;
1045    if(trackaod0->P()<1.){
1046       fPidObjprot->SetTOF(kTRUE);
1047       fPidObjprot->SetTOFdecide(kTRUE);
1048    }
1049   }
1050  }else{
1051    if(track0->P()<1.2){
1052     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1053     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1054    }else{
1055     AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(0));
1056     if(trackaod0->P()<1.){
1057       fPidObjprot->SetTOF(kFALSE);
1058       fPidObjprot->SetTOFdecide(kFALSE);
1059     }
1060     Int_t isProton=fPidObjprot->MakeRawPid(trackaod0,4);
1061     if(isProton>=1) isp0=kTRUE;
1062     if(trackaod0->P()<1.){
1063       fPidObjprot->SetTOF(kTRUE);
1064       fPidObjprot->SetTOFdecide(kTRUE);
1065     }
1066    }
1067  }
1068
1069  if(isTOF2){ //proton
1070   if(track2->P()<2.2) {
1071     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
1072     if(obj->Pt()<3. && track2->P()<1.) fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1073     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1074   }else{
1075     AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1076     if(trackaod2->P()<1.){
1077       fPidObjprot->SetTOF(kFALSE);
1078       fPidObjprot->SetTOFdecide(kFALSE);
1079     }
1080    Int_t isProton=fPidObjprot->MakeRawPid(trackaod2,4);
1081    if(isProton>=1) isp2=kTRUE;
1082    if(trackaod2->P()<1.){
1083       fPidObjprot->SetTOF(kTRUE);
1084       fPidObjprot->SetTOFdecide(kTRUE);
1085     }
1086   }
1087  }else{
1088    if(track2->P()<1.2){
1089     fPidHF->GetPidCombined()->SetDetectorMask(AliPIDResponse::kDetTPC);
1090     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1091    }else{
1092     AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1093      if(trackaod2->P()<1.){
1094       fPidObjprot->SetTOF(kFALSE);
1095       fPidObjprot->SetTOFdecide(kFALSE);
1096     }
1097     Int_t isProton=fPidObjprot->MakeRawPid(trackaod2,4);
1098     if(isProton>=1) isp2=kTRUE;
1099     if(trackaod2->P()<1.){
1100       fPidObjprot->SetTOF(kTRUE);
1101       fPidObjprot->SetTOFdecide(kTRUE);
1102     }
1103    }
1104  }
1105   AliAODTrack *trackaod2=(AliAODTrack*)(obj->GetDaughter(2));
1106   if(fPidObjpion->MakeRawPid(trackaod2,2)>=1)ispi2=kTRUE;
1107   AliAODTrack *trackaod0=(AliAODTrack*)(obj->GetDaughter(2));
1108   if(fPidObjpion->MakeRawPid(trackaod0,2)>=1)ispi0=kTRUE;
1109
1110   if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
1111
1112     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;
1113     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1;
1114     }
1115
1116    if(!isK1 && TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]) isK1=kTRUE;
1117     if(!ispi0 && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) ispi0=kTRUE;
1118     if(!ispi2 && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) ispi2=kTRUE;
1119     if(!isp0 && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton]) isp0=kTRUE;
1120     if(!isp2 && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton]) isp2=kTRUE;
1121     if(isK1 && ispi0 && isp2) okLcpiKp = 1;
1122     if(isK1 && isp0 && ispi2) okLcpKpi = 1;
1123
1124    if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1125     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1126     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1127
1128     return returnvalue;
1129
1130
1131 }
1132 //----------------------------------------------------------
1133 Int_t AliRDHFCutsLctopKpi::IsSelectedPIDStrong(AliAODRecoDecayHF* obj) {
1134
1135
1136     if(!fUsePID || !obj) return 3;
1137     Int_t okLcpKpi=0,okLcpiKp=0;
1138     Int_t returnvalue=0;
1139     Bool_t isPeriodd=fPidHF->GetOnePad();
1140     Bool_t isMC=fPidHF->GetMC();
1141     Bool_t ispion0=kTRUE,ispion2=kTRUE;
1142     Bool_t isproton0=kFALSE,isproton2=kFALSE;
1143     Bool_t iskaon1=kFALSE;
1144     if(isPeriodd) {
1145      fPidObjprot->SetOnePad(kTRUE);
1146      fPidObjpion->SetOnePad(kTRUE);
1147     }
1148     if(isMC) {
1149      fPidObjprot->SetMC(kTRUE);
1150      fPidObjpion->SetMC(kTRUE);
1151     }
1152
1153     for(Int_t i=0;i<3;i++){
1154      AliAODTrack *track=(AliAODTrack*)obj->GetDaughter(i);
1155      if(!track) return 0;
1156      // identify kaon
1157      if(i==1) {
1158       Int_t isKaon=fPidHF->MakeRawPid(track,3);
1159       if(isKaon>=1) {
1160        iskaon1=kTRUE;
1161       if(fPidHF->MakeRawPid(track,2)>=1) iskaon1=kFALSE;
1162       }
1163       if(!iskaon1) return 0;
1164
1165      }else{
1166      //pion or proton
1167
1168      Int_t isProton=fPidObjprot->MakeRawPid(track,4);
1169      if(isProton>=1){
1170       if(fPidHF->MakeRawPid(track,2)>=1) isProton=-1;
1171       if(fPidHF->MakeRawPid(track,3)>=1) isProton=-1;
1172      }
1173
1174      Int_t isPion=fPidObjpion->MakeRawPid(track,2);
1175      if(fPidHF->MakeRawPid(track,3)>=1) isPion=-1;
1176      if(fPidObjprot->MakeRawPid(track,4)>=1) isPion=-1;
1177
1178
1179      if(i==0) {
1180       if(isPion<0) ispion0=kFALSE;
1181       if(isProton>=1) isproton0=kTRUE;
1182
1183      }
1184       if(!ispion0 && !isproton0) return 0;
1185      if(i==2) {
1186       if(isPion<0) ispion2=kFALSE;
1187       if(isProton>=1) isproton2=kTRUE;
1188      }
1189
1190     }
1191    }
1192
1193     if(ispion2 && isproton0 && iskaon1) okLcpKpi=1;
1194     if(ispion0 && isproton2 && iskaon1) okLcpiKp=1;
1195     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1196     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1197     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1198
1199  return returnvalue;
1200 }
1201 //--------------------
1202 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDpPb(AliAODRecoDecayHF* obj) {
1203
1204     if(!fUsePID || !obj) {return 3;}
1205     Int_t okLcpKpi=0,okLcpiKp=0;
1206     Int_t returnvalue=0;
1207
1208     Bool_t isMC=fPidHF->GetMC();
1209
1210
1211     if(isMC) {
1212             fPidObjprot->SetMC(kTRUE);
1213             fPidObjpion->SetMC(kTRUE);
1214     }
1215
1216     AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
1217     AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
1218     AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
1219     if (!track0 || !track1 || !track2) return 0;
1220     Double_t prob0[AliPID::kSPECIES];
1221     Double_t prob1[AliPID::kSPECIES];
1222     Double_t prob2[AliPID::kSPECIES];
1223
1224     fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1225     fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
1226     fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1227
1228
1229     if(fPIDThreshold[AliPID::kPion]>0. && fPIDThreshold[AliPID::kKaon]>0. && fPIDThreshold[AliPID::kProton]>0.){
1230     okLcpiKp=  (prob0[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ])
1231              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
1232              &&(prob2[AliPID::kProton]>fPIDThreshold[AliPID::kProton]);
1233     okLcpKpi=  (prob0[AliPID::kProton]>fPIDThreshold[AliPID::kProton])
1234              &&(prob1[AliPID::kKaon  ]>fPIDThreshold[AliPID::kKaon  ])
1235              &&(prob2[AliPID::kPion  ]>fPIDThreshold[AliPID::kPion  ]);
1236    }else{
1237                     //pion or proton
1238
1239
1240     if(TMath::MaxElement(AliPID::kSPECIES,prob1) == prob1[AliPID::kKaon]){
1241     if(TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kPion]) okLcpKpi = 1;
1242     if(TMath::MaxElement(AliPID::kSPECIES,prob2) == prob2[AliPID::kProton] && TMath::MaxElement(AliPID::kSPECIES,prob0) == prob0[AliPID::kPion]) okLcpiKp = 1;
1243             }
1244            }
1245
1246     if(okLcpKpi) returnvalue=1; //cuts passed as Lc->pKpi
1247     if(okLcpiKp) returnvalue=2; //cuts passed as Lc->piKp
1248     if(okLcpKpi && okLcpiKp) returnvalue=3; //cuts passed as both pKpi and piKp
1249
1250     return returnvalue;
1251 }
1252 //-----------------------
1253 Int_t AliRDHFCutsLctopKpi::IsSelectedCombinedPIDpPb2(AliAODRecoDecayHF* obj) {
1254
1255   Int_t returnvalue =0;
1256   Double_t thresholdK =0.7;
1257   Double_t thresholdPi =0.3; //!
1258   Double_t thresholdPr =0.7;
1259
1260   AliVTrack *track0=dynamic_cast<AliVTrack*>(obj->GetDaughter(0));
1261   AliVTrack *track1=dynamic_cast<AliVTrack*>(obj->GetDaughter(1));
1262   AliVTrack *track2=dynamic_cast<AliVTrack*>(obj->GetDaughter(2));
1263
1264   if (!track0 || !track1 || !track2) return 0;
1265   Double_t prob0[AliPID::kSPECIES];
1266   Double_t prob1[AliPID::kSPECIES];
1267   Double_t prob2[AliPID::kSPECIES];
1268
1269   fPidHF->GetPidCombined()->ComputeProbabilities(track0,fPidHF->GetPidResponse(),prob0);
1270   fPidHF->GetPidCombined()->ComputeProbabilities(track1,fPidHF->GetPidResponse(),prob1);
1271   fPidHF->GetPidCombined()->ComputeProbabilities(track2,fPidHF->GetPidResponse(),prob2);
1272
1273   if(prob1[AliPID::kKaon]>thresholdK){
1274     if(TMath::MaxElement(AliPID::kSPECIES,prob0)>TMath::MaxElement(AliPID::kSPECIES,prob2)){
1275       if(((prob0[AliPID::kPion  ]>prob0[AliPID::kProton  ])&& prob0[AliPID::kPion]>thresholdPi) && prob2[AliPID::kProton]>thresholdPr) returnvalue=2;//piKp
1276       else if(((prob0[AliPID::kProton  ]>prob0[AliPID::kPion  ])&& prob0[AliPID::kProton]>thresholdPr) && prob2[AliPID::kPion]>thresholdPi)returnvalue =1;//pKpi
1277     }
1278
1279     else if(TMath::MaxElement(AliPID::kSPECIES,prob0)<TMath::MaxElement(AliPID::kSPECIES,prob2)){
1280       if(((prob2[AliPID::kPion  ]>prob2[AliPID::kProton  ])&& prob2[AliPID::kPion]>thresholdPi) && prob0[AliPID::kProton]>thresholdPr) returnvalue=1; //pKpi
1281       else if(((prob2[AliPID::kProton  ]>prob2[AliPID::kPion  ])&& prob2[AliPID::kProton]>thresholdPr) && prob0[AliPID::kPion]>thresholdPi)returnvalue =2;      //piKp
1282     }
1283
1284   }
1285   return returnvalue;
1286
1287 }
1288 //------------------------------------------------------
1289 void AliRDHFCutsLctopKpi::SetStandardCutsPPb2013() {
1290
1291  SetName("LctopKpiProdCuts");
1292  SetTitle("Production cuts for Lc analysis");
1293
1294  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
1295  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
1296  esdTrackCuts->SetRequireTPCRefit(kTRUE);
1297  esdTrackCuts->SetMinNClustersTPC(70);
1298  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,
1299                                           AliESDtrackCuts::kAny);
1300  esdTrackCuts->SetRequireITSRefit(kTRUE);
1301  esdTrackCuts->SetMinNClustersITS(4);
1302  esdTrackCuts->SetMinDCAToVertexXY(0.);
1303  esdTrackCuts->SetEtaRange(-0.8,0.8);
1304  esdTrackCuts->SetPtRange(0.3,1.e10);
1305  AddTrackCuts(esdTrackCuts);
1306  delete esdTrackCuts;
1307  esdTrackCuts=NULL;
1308
1309  const Int_t nvars=13;
1310  const Int_t nptbins=9;
1311  Float_t* ptbins;
1312  ptbins=new Float_t[nptbins+1];
1313  ptbins[0]=0.;
1314  ptbins[1]=1.;
1315  ptbins[2]=2.;
1316  ptbins[3]=3.;
1317  ptbins[4]=4.;
1318  ptbins[5]=5.;
1319  ptbins[6]=6.;
1320  ptbins[7]=8.;
1321  ptbins[8]=10.;
1322  ptbins[9]=99999.;
1323
1324  SetGlobalIndex(nvars,nptbins);
1325  SetPtBins(nptbins+1,ptbins);
1326
1327  Float_t** prodcutsval;
1328  prodcutsval=new Float_t*[nvars];
1329  for(Int_t iv=0;iv<nvars;iv++){
1330   prodcutsval[iv]=new Float_t[nptbins];
1331  }
1332
1333  for(Int_t ipt=0;ipt<nptbins;ipt++){
1334     prodcutsval[0][ipt]=0.13;
1335     prodcutsval[1][ipt]=0.4;
1336     prodcutsval[2][ipt]=0.4;
1337     prodcutsval[3][ipt]=0.;
1338     prodcutsval[4][ipt]=0.;
1339     prodcutsval[5][ipt]=0.;
1340     prodcutsval[6][ipt]=0.06;
1341     prodcutsval[7][ipt]=0.;
1342     prodcutsval[8][ipt]=0.;
1343     prodcutsval[9][ipt]=0.;
1344     prodcutsval[10][ipt]=0.;
1345     prodcutsval[11][ipt]=0.05;
1346     prodcutsval[12][ipt]=0.4;
1347  }
1348  SetCuts(nvars,nptbins,prodcutsval);
1349
1350  AliAODPidHF* pidObjK=new AliAODPidHF();
1351  Double_t sigmasK[5]={3.,1.,1.,3.,2.};
1352  pidObjK->SetSigma(sigmasK);
1353  pidObjK->SetAsym(kTRUE);
1354  pidObjK->SetMatch(1);
1355  pidObjK->SetTPC(kTRUE);
1356  pidObjK->SetTOF(kTRUE);
1357  pidObjK->SetITS(kTRUE);
1358  Double_t plimK[2]={0.5,0.8};
1359  pidObjK->SetPLimit(plimK,2);
1360  pidObjK->SetTOFdecide(kTRUE);
1361  SetPidHF(pidObjK);
1362
1363  AliAODPidHF* pidObjpi=new AliAODPidHF();
1364  pidObjpi->SetTPC(kTRUE);
1365  Double_t sigmaspi[5]={3.,0.,0.,0.,0.};
1366  pidObjpi->SetSigma(sigmaspi);
1367  pidObjpi->SetTOFdecide(kTRUE);
1368  SetPidpion(pidObjpi);
1369
1370  AliAODPidHF* pidObjp=new AliAODPidHF();
1371  Double_t sigmasp[5]={3.,1.,1.,3.,2.};
1372  pidObjp->SetSigma(sigmasp);
1373  pidObjp->SetAsym(kTRUE);
1374  pidObjp->SetMatch(1);
1375  pidObjp->SetTPC(kTRUE);
1376  pidObjp->SetTOF(kTRUE);
1377  pidObjp->SetITS(kTRUE);
1378  Double_t plimp[2]={1.,2.};
1379  pidObjp->SetPLimit(plimp,2);
1380  pidObjp->SetTOFdecide(kTRUE);
1381
1382  SetUsePID(kTRUE);
1383
1384  PrintAll();
1385
1386  for(Int_t iiv=0;iiv<nvars;iiv++){
1387   delete [] prodcutsval[iiv];
1388  }
1389  delete [] prodcutsval;
1390  prodcutsval=NULL;
1391  delete [] ptbins;
1392  ptbins=NULL;
1393
1394  delete pidObjK;
1395  pidObjK=NULL;
1396  delete pidObjpi;
1397  pidObjpi=NULL;
1398  delete pidObjp;
1399  pidObjp=NULL;
1400
1401  return;
1402 }
1403