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