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