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