]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliRDHFCutsD0toKpi.cxx
Some fixes
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliRDHFCutsD0toKpi.cxx
index 549a0acedf398f0cbca68d9b86103631090a8dfd..5729b603120cec9bf8e1ebc394435d881e1e0892 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 /////////////////////////////////////////////////////////////
 //
 // Class for cuts on AOD reconstructed D0->Kpi
 #include "AliAODPid.h"
 #include "AliTPCPIDResponse.h"
 #include "AliAODVertex.h"
+#include "AliKFVertex.h"
+#include "AliKFParticle.h"
 
 ClassImp(AliRDHFCutsD0toKpi)
 
 //--------------------------------------------------------------------------
 AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const char* name) : 
 AliRDHFCuts(name),
-fUseSpecialCuts(kFALSE),
+fUseSpecialCuts(kTRUE),
 fLowPt(kTRUE),
-fDefaultPID(kTRUE)
+fDefaultPID(kFALSE),
+fUseKF(kFALSE)
 {
   //
   // Default Constructor
   //
-  Int_t nvars=9;
+  Int_t nvars=11;
   SetNVars(nvars);
-  TString varNames[9]={"inv. mass [GeV]",   
-                      "dca [cm]",
-                      "cosThetaStar", 
-                      "pTK [GeV/c]",
-                      "pTPi [GeV/c]",
-                      "d0K [cm]",
-                      "d0Pi [cm]",
-                      "d0d0 [cm^2]",
-                      "cosThetaPoint"};
-  Bool_t isUpperCut[9]={kTRUE,
-                       kTRUE,
-                       kTRUE,
-                       kFALSE,
-                       kFALSE,
-                       kTRUE,
-                       kTRUE,
-                       kTRUE,
-                       kFALSE};
+  TString varNames[11]={"inv. mass [GeV]",   
+                       "dca [cm]",
+                       "cosThetaStar", 
+                       "pTK [GeV/c]",
+                       "pTPi [GeV/c]",
+                       "d0K [cm]",
+                       "d0Pi [cm]",
+                       "d0d0 [cm^2]",
+                       "cosThetaPoint",
+                       "|cosThetaPointXY|", 
+                       "NormDecayLenghtXY"};
+  Bool_t isUpperCut[11]={kTRUE,
+                        kTRUE,
+                        kTRUE,
+                        kFALSE,
+                        kFALSE,
+                        kTRUE,
+                        kTRUE,
+                        kTRUE,
+                        kFALSE,
+                        kFALSE, 
+                        kFALSE};
   SetVarNames(nvars,varNames,isUpperCut);
-  Bool_t forOpt[9]={kFALSE,
-                   kTRUE,
-                   kTRUE,
-                   kFALSE,
-                   kFALSE,
-                   kFALSE,
-                   kFALSE,
-                   kTRUE,
-                   kTRUE};
+  Bool_t forOpt[11]={kFALSE,
+                    kTRUE,
+                    kTRUE,
+                    kFALSE,
+                    kFALSE,
+                    kFALSE,
+                    kFALSE,
+                    kTRUE,
+                    kTRUE,
+                    kFALSE,
+                    kFALSE};
   SetVarsForOpt(4,forOpt);
   Float_t limits[2]={0,999999999.};
   SetPtBins(2,limits);
@@ -83,7 +94,8 @@ AliRDHFCutsD0toKpi::AliRDHFCutsD0toKpi(const AliRDHFCutsD0toKpi &source) :
   AliRDHFCuts(source),   
   fUseSpecialCuts(source.fUseSpecialCuts),
   fLowPt(source.fLowPt),
-  fDefaultPID(source.fDefaultPID)
+  fDefaultPID(source.fDefaultPID),
+  fUseKF(source.fUseKF)
 {
   //
   // Copy constructor
@@ -108,7 +120,7 @@ AliRDHFCutsD0toKpi &AliRDHFCutsD0toKpi::operator=(const AliRDHFCutsD0toKpi &sour
 
 
 //---------------------------------------------------------------------------
-void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters) {
+void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int_t nvars,Int_t *pdgdaughters,AliAODEvent *aod) {
   // 
   // Fills in vars the values of the variables 
   //
@@ -120,6 +132,18 @@ void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int
 
   AliAODRecoDecayHF2Prong *dd = (AliAODRecoDecayHF2Prong*)d;
  
+  //recalculate vertex w/o daughters
+  Bool_t cleanvtx=kFALSE;
+  AliAODVertex *origownvtx=0x0;
+  if(fRemoveDaughtersFromPrimary) {
+    if(dd->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*dd->GetOwnPrimaryVtx());
+    cleanvtx=kTRUE;
+    if(!RecalcOwnPrimaryVtx(dd,aod)) {
+      CleanOwnPrimaryVtx(dd,aod,origownvtx);
+      cleanvtx=kFALSE;
+    }
+  }
+
   Int_t iter=-1;
   if(fVarsForOpt[0]){
     iter++;
@@ -186,6 +210,18 @@ void AliRDHFCutsD0toKpi::GetCutVarsForOpt(AliAODRecoDecayHF *d,Float_t *vars,Int
     vars[iter]=dd->CosPointingAngle();
   }
   
+  if(fVarsForOpt[9]){
+               iter++;
+               vars[iter]=TMath::Abs(dd->CosPointingAngleXY());
+       }
+  
+   if(fVarsForOpt[10]){
+               iter++;
+          vars[iter]=(dd->NormalizedDecayLengthXY()*(dd->P()/dd->Pt()));
+       }
+
+   if(cleanvtx)CleanOwnPrimaryVtx(dd,aod,origownvtx);
+
   return;
 }
 //---------------------------------------------------------------------------
@@ -193,6 +229,10 @@ Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEve
   //
   // Apply selection
   //
+
+  fIsSelectedCuts=0;
+  fIsSelectedPID=0;
 
   if(!fCutsRD){
     cout<<"Cut matrice not inizialized. Exit..."<<endl;
@@ -206,132 +246,386 @@ Int_t AliRDHFCutsD0toKpi::IsSelected(TObject* obj,Int_t selectionLevel,AliAODEve
     return 0;
   }
 
-  // selection on daughter tracks 
-  if(selectionLevel==AliRDHFCuts::kAll || 
-     selectionLevel==AliRDHFCuts::kTracks) {
-    if(!AreDaughtersSelected(d)) return 0;
-  }
-
-
+  Double_t ptD=d->Pt();
+  if(ptD<fMinPtCand) return 0;
+  if(ptD>fMaxPtCand) return 0;
  
   // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
   Int_t returnvaluePID=3;
   Int_t returnvalueCuts=3;
 
-  // selection on PID 
-  if(selectionLevel==AliRDHFCuts::kAll || 
-     selectionLevel==AliRDHFCuts::kCandidate ||
-     selectionLevel==AliRDHFCuts::kPID) {
-    returnvaluePID = IsSelectedPID(d);
-  }
-
-
-
   // selection on candidate
   if(selectionLevel==AliRDHFCuts::kAll || 
      selectionLevel==AliRDHFCuts::kCandidate) {
 
-    //recalculate vertex w/o daughters
-    AliAODVertex *origownvtx=0x0;
-    AliAODVertex *recvtx=0x0;
-  
-    if(fRemoveDaughtersFromPrimary) {
-      if(!aod) {
-       AliError("Can not remove daughters from vertex without AOD event");
-       return 0;
-      }   
-      if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
-      recvtx=d->RemoveDaughtersFromPrimaryVtx(aod);
-      if(!recvtx){
-       AliDebug(2,"Removal of daughter tracks failed");
-       //recvtx=d->GetPrimaryVtx();
-       if(origownvtx){
-         delete origownvtx;
-         origownvtx=NULL;
+    if(!fUseKF) {
+
+      //recalculate vertex w/o daughters
+      AliAODVertex *origownvtx=0x0;
+      if(fRemoveDaughtersFromPrimary && !fUseMCVertex) {
+       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+       if(!RecalcOwnPrimaryVtx(d,aod)) { 
+         CleanOwnPrimaryVtx(d,aod,origownvtx);
+         return 0;
        }
-       return 0;
       }
-      //set recalculed primary vertex
-      d->SetOwnPrimaryVtx(recvtx);
-      delete recvtx; recvtx=NULL;
-    }
 
-    
-    Double_t pt=d->Pt();
-   
-    Int_t okD0=0,okD0bar=0;
-    Int_t ptbin=PtBin(pt);
-    if (ptbin==-1) {
-      if(origownvtx){
-       d->SetOwnPrimaryVtx(origownvtx);
-       delete origownvtx;
-       origownvtx=NULL;
+      if(fUseMCVertex) {
+       if(d->GetOwnPrimaryVtx()) origownvtx=new AliAODVertex(*d->GetOwnPrimaryVtx());
+       if(!SetMCPrimaryVtx(d,aod)) {
+         CleanOwnPrimaryVtx(d,aod,origownvtx);
+         return 0;
+       }
+      }
+      
+      Double_t pt=d->Pt();
+      
+      Int_t okD0=0,okD0bar=0;
+      
+      Int_t ptbin=PtBin(pt);
+      if (ptbin==-1) {
+       CleanOwnPrimaryVtx(d,aod,origownvtx);
+       return 0;
       }
-      else d->UnsetOwnPrimaryVtx();
-      return 0;
-    }
-    Double_t mD0,mD0bar,ctsD0,ctsD0bar;
-    okD0=1; okD0bar=1;
 
-    Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+      Double_t mD0,mD0bar,ctsD0,ctsD0bar;
+      okD0=1; okD0bar=1;
+      
+      Double_t mD0PDG = TDatabasePDG::Instance()->GetParticle(421)->Mass();
 
-    if(d->PtProng(1) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(0) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
-    if(d->PtProng(0) < fCutsRD[GetGlobalIndex(3,ptbin)] || d->PtProng(1) < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
-    if(!okD0 && !okD0bar) returnvalueCuts=0;
-  
-    
-    if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
-       TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
-    if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
-       TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
-    if(!okD0 && !okD0bar)  returnvalueCuts=0;
-    
-    if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)])  returnvalueCuts=0;
-    
-    d->InvMassD0(mD0,mD0bar);
-    if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
-    if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)])  okD0bar = 0;
-    if(!okD0 && !okD0bar)  returnvalueCuts=0;
-    
-    d->CosThetaStarD0(ctsD0,ctsD0bar);
-    if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; 
-    if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
-    if(!okD0 && !okD0bar)   returnvalueCuts=0;
+      d->InvMassD0(mD0,mD0bar);
+      if(TMath::Abs(mD0-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
+      if(TMath::Abs(mD0bar-mD0PDG) > fCutsRD[GetGlobalIndex(0,ptbin)])  okD0bar = 0;
+      if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+      if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
     
-    if(d->Prodd0d0() > fCutsRD[GetGlobalIndex(7,ptbin)])  returnvalueCuts=0;
+      
+      if(d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
+      if(d->Pt2Prong(0) < fCutsRD[GetGlobalIndex(3,ptbin)]*fCutsRD[GetGlobalIndex(3,ptbin)] || d->Pt2Prong(1) < fCutsRD[GetGlobalIndex(4,ptbin)]*fCutsRD[GetGlobalIndex(4,ptbin)]) okD0bar = 0;
+      if(!okD0 && !okD0bar) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+      
+      
+      if(TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
+        TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
+      if(TMath::Abs(d->Getd0Prong(0)) > fCutsRD[GetGlobalIndex(6,ptbin)] ||
+        TMath::Abs(d->Getd0Prong(1)) > fCutsRD[GetGlobalIndex(5,ptbin)]) okD0bar = 0;
+      if(!okD0 && !okD0bar)  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+      
+      if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+      
     
-    if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)])  returnvalueCuts=0;
+      d->CosThetaStarD0(ctsD0,ctsD0bar);
+      if(TMath::Abs(ctsD0) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; 
+      if(TMath::Abs(ctsD0bar) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0bar = 0;
+      if(!okD0 && !okD0bar)   {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
     
-    if (returnvalueCuts!=0) {
-      if (okD0) returnvalueCuts=1; //cuts passed as D0
-      if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
-      if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
-    }
+      if(d->CosPointingAngle() < fCutsRD[GetGlobalIndex(8,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+      
+      if(TMath::Abs(d->CosPointingAngleXY()) < fCutsRD[GetGlobalIndex(9,ptbin)])  {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+       
+      Double_t normalDecayLengXY=(d->NormalizedDecayLengthXY()*(d->P()/d->Pt()));
+      if (normalDecayLengXY < fCutsRD[GetGlobalIndex(10, ptbin)]) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+      
+      if (returnvalueCuts!=0) {
+       if (okD0) returnvalueCuts=1; //cuts passed as D0
+       if (okD0bar) returnvalueCuts=2; //cuts passed as D0bar
+       if (okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
+      }
+
+      // call special cuts
+      Int_t special=1;
+      if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
+      if(!special) {CleanOwnPrimaryVtx(d,aod,origownvtx); return 0;}
+
+      // unset recalculated primary vertex when not needed any more
+      CleanOwnPrimaryVtx(d,aod,origownvtx);
 
-    // call special cuts
-    Int_t special=1;
-    if(fUseSpecialCuts) special=IsSelectedSpecialCuts(d);
-    if(!special) returnvalueCuts=0;
-
-    // unset recalculated primary vertex when not needed any more
-    if(origownvtx) {
-      d->SetOwnPrimaryVtx(origownvtx);
-      delete origownvtx;
-      origownvtx=NULL;
-    } else if(fRemoveDaughtersFromPrimary) {
-      d->UnsetOwnPrimaryVtx();
-      AliDebug(3,"delete new vertex\n");
+    } else {
+      // go to selection with Kalman vertexing, if requested
+      returnvalueCuts = IsSelectedKF(d,aod);
     }
 
+    fIsSelectedCuts=returnvalueCuts;
+    if(!returnvalueCuts) return 0;
   }
 
 
+  // selection on PID 
+  if(selectionLevel==AliRDHFCuts::kAll || 
+     selectionLevel==AliRDHFCuts::kCandidate ||
+     selectionLevel==AliRDHFCuts::kPID) {
+    returnvaluePID = IsSelectedPID(d);
+    fIsSelectedPID=returnvaluePID;
+    if(!returnvaluePID) return 0;
+  }
+
+  Int_t returnvalueComb=CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
+
+  if(!returnvalueComb) return 0;
+
+  // selection on daughter tracks 
+  if(selectionLevel==AliRDHFCuts::kAll || 
+     selectionLevel==AliRDHFCuts::kTracks) {
+    if(!AreDaughtersSelected(d)) return 0;
+  }
   //  cout<<"Pid = "<<returnvaluePID<<endl;
-  return CombineSelectionLevels(3,returnvalueCuts,returnvaluePID);
+  return returnvalueComb;
 }
+
+//------------------------------------------------------------------------------------------
+Int_t AliRDHFCutsD0toKpi::IsSelectedKF(AliAODRecoDecayHF2Prong *d,
+                                      AliAODEvent* aod) const {
+  //
+  // Apply selection using KF-vertexing
+  //
+       
+  AliAODTrack *track0 = (AliAODTrack*)d->GetDaughter(0);
+  AliAODTrack *track1 = (AliAODTrack*)d->GetDaughter(1); 
+   
+  if(!track0 || !track1) {
+    cout<<"one or two D0 daughters missing!"<<endl;
+    return 0;
+  }
+
+  // returnvalue: 0 not sel, 1 only D0, 2 only D0bar, 3 both
+  Int_t returnvalueCuts=3;
+        
+  // candidate selection with AliKF
+  AliKFParticle::SetField(aod->GetMagneticField()); // set the magnetic field
+  
+  Int_t okD0=0,okD0bar=0;
+  okD0=1; okD0bar=1;
+  
+  // convert tracks into AliKFParticles
+  
+  AliKFParticle negPiKF(*track1,-211); // neg pion kandidate
+  AliKFParticle negKKF(*track1,-321); // neg kaon kandidate
+  AliKFParticle posPiKF(*track0,211); // pos pion kandidate
+  AliKFParticle posKKF(*track0,321); // pos kaon kandidate
+  
+  // build D0 candidates
+  
+  AliKFParticle d0c(negKKF,posPiKF); // D0 candidate
+  AliKFParticle ad0c(posKKF,negPiKF); // D0bar candidate
+  
+  // create kf primary vertices
+  
+  AliAODVertex *vtx1 = aod->GetPrimaryVertex();
+  AliKFVertex primVtx(*vtx1); 
+  AliKFVertex aprimVtx(*vtx1);
+  
+  if(primVtx.GetNContributors()<=0) okD0 = 0;
+  if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
+  if(!okD0 && !okD0bar) returnvalueCuts=0;
+       
+  // calculate mass
+       
+  Double_t d0mass = d0c.GetMass();
+  Double_t ad0mass = ad0c.GetMass();
+       
+  // calculate P of D0 and D0bar
+  Double_t d0P = d0c.GetP();
+  Double_t d0Px = d0c.GetPx();
+  Double_t d0Py = d0c.GetPy();
+  Double_t d0Pz = d0c.GetPz();
+  Double_t ad0P = ad0c.GetP(); 
+  Double_t ad0Px = ad0c.GetPx();
+  Double_t ad0Py = ad0c.GetPy();
+  Double_t ad0Pz = ad0c.GetPz();
+  
+  //calculate Pt of D0 and D0bar
+       
+  Double_t pt=d0c.GetPt(); 
+  Double_t apt=ad0c.GetPt();
+       
+  // remove D0 daughters from primary vertices (if used in vertex fit) and add D0-candidates
+  
+  if(track0->GetUsedForPrimVtxFit()) {
+    primVtx -= posPiKF; 
+    aprimVtx -= posKKF;
+  }
+  
+  if(track1->GetUsedForPrimVtxFit()) { 
+    primVtx -= negKKF; 
+    aprimVtx -= negPiKF;
+  }
+  
+  primVtx += d0c;
+  aprimVtx += ad0c;
+  
+  if(primVtx.GetNContributors()<=0) okD0 = 0;
+  if(aprimVtx.GetNContributors()<=0) okD0bar = 0;
+  if(!okD0 && !okD0bar) returnvalueCuts=0;
+  
+  //calculate cut variables
+  
+  // calculate impact params of daughters w.r.t recalculated vertices
+  
+  Double_t impactPi = posPiKF.GetDistanceFromVertexXY(primVtx);
+  Double_t aimpactPi = negPiKF.GetDistanceFromVertexXY(aprimVtx);
+  Double_t impactKa = negKKF.GetDistanceFromVertexXY(primVtx);
+  Double_t aimpactKa = posKKF.GetDistanceFromVertexXY(aprimVtx);
+       
+  // calculate Product of Impact Params
+       
+  Double_t prodParam = impactPi*impactKa;
+  Double_t aprodParam = aimpactPi*aimpactKa;
+       
+  // calculate cosine of pointing angles
+       
+  TVector3 mom(d0c.GetPx(),d0c.GetPy(),d0c.GetPz());
+  TVector3 fline(d0c.GetX()-primVtx.GetX(),
+                d0c.GetY()-primVtx.GetY(),
+                d0c.GetZ()-primVtx.GetZ());
+  Double_t pta = mom.Angle(fline);
+  Double_t cosP = TMath::Cos(pta); // cosine of pta for D0 candidate
+  
+  TVector3 amom(ad0c.GetPx(),ad0c.GetPy(),ad0c.GetPz());
+  TVector3 afline(ad0c.GetX()-aprimVtx.GetX(),
+                 ad0c.GetY()-aprimVtx.GetY(),
+                 ad0c.GetZ()-aprimVtx.GetZ());
+  Double_t apta = amom.Angle(afline);
+  Double_t acosP = TMath::Cos(apta); // cosine of pta for D0bar candidate
+  
+  // calculate P of Pions at Decay Position of D0 and D0bar candidates
+  negKKF.TransportToParticle(d0c);
+  posPiKF.TransportToParticle(d0c);
+  posKKF.TransportToParticle(ad0c);
+  negPiKF.TransportToParticle(ad0c);
+  
+  Double_t pxPi =  posPiKF.GetPx();
+  Double_t pyPi =  posPiKF.GetPy();
+  Double_t pzPi =  posPiKF.GetPz();
+  Double_t ptPi =  posPiKF.GetPt();
+  
+  Double_t apxPi =  negPiKF.GetPx();
+  Double_t apyPi =  negPiKF.GetPy();
+  Double_t apzPi =  negPiKF.GetPz();
+  Double_t aptPi =  negPiKF.GetPt();
+  
+  // calculate Pt of Kaons at Decay Position of D0 and D0bar candidates
+  
+  Double_t ptK =  negKKF.GetPt();
+  Double_t aptK =  posKKF.GetPt();
+       
+  //calculate cos(thetastar)
+  Double_t massvtx = TDatabasePDG::Instance()->GetParticle(421)->Mass();
+  Double_t massp[2];
+  massp[0] = TDatabasePDG::Instance()->GetParticle(321)->Mass();
+  massp[1] = TDatabasePDG::Instance()->GetParticle(211)->Mass();
+  Double_t pStar = TMath::Sqrt(TMath::Power(massvtx*massvtx-massp[0]*massp[0]-massp[1]*massp[1],2.)
+                              -4.*massp[0]*massp[0]*massp[1]*massp[1])/(2.*massvtx);
+  
+  // cos(thetastar) for D0 and Pion
+  
+  Double_t d0E = TMath::Sqrt(massvtx*massvtx + d0P*d0P);
+  Double_t beta = d0P/d0E;
+  Double_t gamma = d0E/massvtx;
+  TVector3 momPi(pxPi,pyPi,pzPi);
+  TVector3 momTot(d0Px,d0Py,d0Pz);
+  Double_t q1 = momPi.Dot(momTot)/momTot.Mag();
+  Double_t cts = (q1/gamma-beta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar;     
+       
+  // cos(thetastar) for D0bar and Pion 
+  
+  Double_t ad0E = TMath::Sqrt(massvtx*massvtx + ad0P*ad0P);
+  Double_t abeta = ad0P/ad0E;
+  Double_t agamma = ad0E/massvtx;
+  TVector3 amomPi(apxPi,apyPi,apzPi);
+  TVector3 amomTot(ad0Px,ad0Py,ad0Pz);
+  Double_t aq1 = amomPi.Dot(amomTot)/amomTot.Mag();
+  Double_t acts = (aq1/agamma-abeta*TMath::Sqrt(pStar*pStar+massp[1]*massp[1]))/pStar; 
+  
+  // calculate reduced Chi2 for the full D0 fit
+  d0c.SetProductionVertex(primVtx);
+  ad0c.SetProductionVertex(aprimVtx);
+  negKKF.SetProductionVertex(d0c);
+  posPiKF.SetProductionVertex(d0c);
+  posKKF.SetProductionVertex(ad0c);
+  negPiKF.SetProductionVertex(ad0c);
+  d0c.TransportToProductionVertex();
+  ad0c.TransportToProductionVertex();
+       
+  // calculate the decay length
+  Double_t decayLengthD0 = d0c.GetDecayLength();
+  Double_t adecayLengthD0 = ad0c.GetDecayLength();
+  
+  Double_t chi2D0 = 50.;
+  if(d0c.GetNDF() > 0 && d0c.GetChi2() >= 0) {
+    chi2D0 = d0c.GetChi2()/d0c.GetNDF();
+  }
+  
+  Double_t achi2D0 = 50.;
+  if(ad0c.GetNDF() > 0 && ad0c.GetChi2() >= 0) {
+    achi2D0 = ad0c.GetChi2()/ad0c.GetNDF();
+  }
+       
+  // Get the Pt-bins
+  Int_t ptbin=PtBin(pt);
+  Int_t aptbin=PtBin(apt);
+
+  if(ptbin < 0) okD0 = 0;
+  if(aptbin < 0) okD0bar = 0;
+  if(!okD0 && !okD0bar) returnvalueCuts=0;
+  
+  if(ptK < fCutsRD[GetGlobalIndex(3,ptbin)] || ptPi < fCutsRD[GetGlobalIndex(4,ptbin)]) okD0 = 0;
+  if(aptK < fCutsRD[GetGlobalIndex(3,aptbin)] || aptPi < fCutsRD[GetGlobalIndex(4,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar) returnvalueCuts=0;
+  
+  
+  if(TMath::Abs(impactKa) > fCutsRD[GetGlobalIndex(5,ptbin)] || 
+     TMath::Abs(impactPi) > fCutsRD[GetGlobalIndex(6,ptbin)]) okD0 = 0;
+  
+  if(TMath::Abs(aimpactKa) > fCutsRD[GetGlobalIndex(5,aptbin)] ||
+     TMath::Abs(aimpactPi) > fCutsRD[GetGlobalIndex(6,aptbin)]) okD0bar = 0;
+  
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+  
+  // for the moment via the standard method due to bug in AliKF
+  if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,ptbin)]) okD0 = 0;
+  if(d->GetDCA() > fCutsRD[GetGlobalIndex(1,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar) returnvalueCuts=0;
+    
+    
+  if(TMath::Abs(d0mass-massvtx) > fCutsRD[GetGlobalIndex(0,ptbin)]) okD0 = 0;
+  if(TMath::Abs(ad0mass-massvtx) > fCutsRD[GetGlobalIndex(0,aptbin)])  okD0bar = 0;
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+  
+  
+  if(TMath::Abs(cts) > fCutsRD[GetGlobalIndex(2,ptbin)]) okD0 = 0; 
+  if(TMath::Abs(acts) > fCutsRD[GetGlobalIndex(2,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar)   returnvalueCuts=0;
+  
+  if(prodParam  > fCutsRD[GetGlobalIndex(7,ptbin)]) okD0 = 0;
+  if(aprodParam > fCutsRD[GetGlobalIndex(7,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+    
+  if(cosP  < fCutsRD[GetGlobalIndex(8,ptbin)]) okD0 = 0; 
+  if(acosP < fCutsRD[GetGlobalIndex(8,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+       
+  if(chi2D0  > fCutsRD[GetGlobalIndex(10,ptbin)]) okD0 = 0; 
+  if(achi2D0 > fCutsRD[GetGlobalIndex(10,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+       
+  if(decayLengthD0  < fCutsRD[GetGlobalIndex(9,ptbin)]) okD0 = 0; 
+  if(adecayLengthD0 < fCutsRD[GetGlobalIndex(9,aptbin)]) okD0bar = 0;
+  if(!okD0 && !okD0bar)  returnvalueCuts=0;
+    
+  if(returnvalueCuts!=0) {
+    if(okD0) returnvalueCuts=1; //cuts passed as D0
+    if(okD0bar) returnvalueCuts=2; //cuts passed as D0bar
+    if(okD0 && okD0bar) returnvalueCuts=3; //cuts passed as D0 and D0bar
+  }
+
+  return returnvalueCuts;  
+}
+
 //---------------------------------------------------------------------------
+
 Bool_t AliRDHFCutsD0toKpi::IsInFiducialAcceptance(Double_t pt, Double_t y) const
 {
   //
@@ -387,10 +681,11 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
   for(Int_t daught=0;daught<2;daught++){
     //Loop con prongs
     AliAODTrack *aodtrack=(AliAODTrack*)d->GetDaughter(daught);
+    if(fPidHF->IsTOFPiKexcluded(aodtrack,5.)) return 0; 
     
     if(!(fPidHF->CheckStatus(aodtrack,"TPC")) && !(fPidHF->CheckStatus(aodtrack,"TOF"))) {
-    checkPIDInfo[daught]=kFALSE; 
-    continue;
+      checkPIDInfo[daught]=kFALSE; 
+      continue;
     }
 
     // identify kaon
@@ -408,26 +703,27 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedPID(AliAODRecoDecayHF* d)
      }
 
 
-   if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
-    isD0D0barPID[0]=0;
-    isD0D0barPID[1]=0;
-   }
-   else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
-    if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
-    else isD0D0barPID[0]=0;// if K+ D0 excluded
-   }
-   else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
-    isD0D0barPID[0]=0;
-    isD0D0barPID[1]=0;
-   }
-   else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){ 
-   if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
-   else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
-        }
-   else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
-    if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
-    else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
-   }
+    if(combinedPID[daught][0]<=-1&&combinedPID[daught][1]<=-1){ // if not a K- and not a pi- both D0 and D0bar excluded
+      isD0D0barPID[0]=0;
+      isD0D0barPID[1]=0;
+    }
+    else if(combinedPID[daught][0]==2&&combinedPID[daught][1]>=1){
+      if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;//if K- D0bar excluded
+      else isD0D0barPID[0]=0;// if K+ D0 excluded
+    }
+    /*    else if(combinedPID[daught][0]==1&&combinedPID[daught][1]>=1){
+         isD0D0barPID[0]=0;
+         isD0D0barPID[1]=0;
+         }
+    */
+    else if(combinedPID[daught][0]>=1||combinedPID[daught][1]<=-1){ 
+      if(aodtrack->Charge()==-1)isD0D0barPID[1]=0;// not a D0bar if K- or if pi- excluded
+      else isD0D0barPID[0]=0;//  not a D0 if K+ or if pi+ excluded
+    }
+    else if(combinedPID[daught][0]<=-1||combinedPID[daught][1]>=1){
+      if(aodtrack->Charge()==-1)isD0D0barPID[0]=0;// not a D0 if pi- or if K- excluded
+      else isD0D0barPID[1]=0;// not a D0bar if pi+ or if K+ excluded
+    }
 
     if(fLowPt && d->Pt()<2.){
      Double_t sigmaTPC[3]={3.,2.,0.};
@@ -622,6 +918,7 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedPIDdefault(AliAODRecoDecayHF* d)
      else{
        tofSig=pid->GetTOFsignal(); 
        pid->GetIntegratedTimes(times);
+       if((tofSig-times[3])>5.*160.)return 0;// PROTON REJECTION
        if(TMath::Abs(tofSig-times[3])>3.*160.){
         isKaonPionTOF[daught][0]=-1;
        }
@@ -774,10 +1071,302 @@ Int_t AliRDHFCutsD0toKpi::IsSelectedSpecialCuts(AliAODRecoDecayHF *d) const
   for(Int_t i=0;i<2/*prongs*/;i++){
     if(TMath::Abs(d->Normalizedd0Prong(i))<normd0Cut) return 0; //normd0Cut not passed
   }
-  if(d->DecayLength()<decLengthCut)  return 0; //decLengthCut not passed
-  if(d->NormalizedDecayLength()<normDecLengthCut)  return 0; //decLengthCut not passed
-    
-
+  if(d->DecayLength2()<decLengthCut*decLengthCut)  return 0; //decLengthCut not passed
+  if(d->NormalizedDecayLength2()<normDecLengthCut*normDecLengthCut)  return 0; //decLengthCut not passed
+       
   return returnvalue;
 }
 
+//----------------------------------------------
+void AliRDHFCutsD0toKpi::SetUseKF(Bool_t useKF)
+{
+  //switch on candidate selection via AliKFparticle
+  if(!useKF) return;
+  if(useKF){
+    fUseKF=useKF;
+    Int_t nvarsKF=11;
+    SetNVars(nvarsKF);
+    TString varNamesKF[11]={"inv. mass [GeV]",   
+                           "dca [cm]",
+                           "cosThetaStar", 
+                           "pTK [GeV/c]",
+                           "pTPi [GeV/c]",
+                           "d0K [cm]",
+                           "d0Pi [cm]",
+                           "d0d0 [cm^2]",
+                           "cosThetaPoint"
+                           "DecayLength[cm]",
+                           "RedChi2"};
+    Bool_t isUpperCutKF[11]={kTRUE,
+                            kTRUE,
+                            kTRUE,
+                            kFALSE,
+                            kFALSE,
+                            kTRUE,
+                            kTRUE,
+                            kTRUE,
+                            kFALSE,
+                            kFALSE,
+                            kTRUE};
+    SetVarNames(nvarsKF,varNamesKF,isUpperCutKF);
+    SetGlobalIndex();
+    Bool_t forOpt[11]={kFALSE,
+                   kTRUE,
+                   kTRUE,
+                   kFALSE,
+                   kFALSE,
+                   kFALSE,
+                   kFALSE,
+                   kTRUE,
+                   kTRUE,
+                   kFALSE,
+                   kFALSE};
+    SetVarsForOpt(4,forOpt);
+  }
+  return;
+}
+
+
+void AliRDHFCutsD0toKpi::SetStandardCutsPP2010() {
+  //
+  //STANDARD CUTS USED FOR 2010 pp analysis 
+  //dca cut will be enlarged soon to 400 micron
+  //
+  
+  SetName("D0toKpiCutsStandard");
+  SetTitle("Standard Cuts for D0 analysis");
+  
+  // PILE UP REJECTION
+  SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+
+  // EVENT CUTS
+  SetMinVtxContr(1);
+
+  
+  // TRACKS ON SINGLE TRACKS
+  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+  esdTrackCuts->SetRequireTPCRefit(kTRUE);
+  esdTrackCuts->SetRequireITSRefit(kTRUE);
+  //  esdTrackCuts->SetMinNClustersITS(4);
+  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+  esdTrackCuts->SetMinDCAToVertexXY(0.);
+  esdTrackCuts->SetEtaRange(-0.8,0.8);
+  esdTrackCuts->SetPtRange(0.3,1.e10);
+  
+  AddTrackCuts(esdTrackCuts);
+  
+  const Int_t nptbins =13;
+  const Double_t ptmax = 9999.;
+  const Int_t nvars=11;
+  Float_t ptbins[nptbins+1];
+  ptbins[0]=0.;
+  ptbins[1]=0.5;       
+  ptbins[2]=1.;
+  ptbins[3]=2.;
+  ptbins[4]=3.;
+  ptbins[5]=4.;
+  ptbins[6]=5.;
+  ptbins[7]=6.;
+  ptbins[8]=8.;
+  ptbins[9]=12.;
+  ptbins[10]=16.;
+  ptbins[11]=20.;
+  ptbins[12]=24.;
+  ptbins[13]=ptmax;
+
+  SetGlobalIndex(nvars,nptbins);
+  SetPtBins(nptbins+1,ptbins);
+  
+  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73,0.,0.},/* pt<0.5*/
+                                                 {0.400,300.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-35000.*1E-8,0.73,0.,0.},/* 0.5<pt<1*/
+                                                 {0.400,200.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-25000.*1E-8,0.75,0.,0.},/* 1<pt<2 */
+                                                 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-15000.*1E-8,0.8,0.,0.},/* 2<pt<3 */
+                                                 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 3<pt<4 */
+                                                 {0.400,200.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 4<pt<5 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 5<pt<6 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-8000.*1E-8,0.85,0.,0.},/* 6<pt<8 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.,0.},/* 8<pt<12 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-0.*1E-8,0.85,0.,0.},/* 12<pt<16 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85,0.,0.},/* 16<pt<20 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85,0.,0.},/* 20<pt<24 */
+                                                 {0.400,150.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,0.*1E-8,0.85,0.,0.}};/* pt>24 */
+  
+  
+  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+  
+  for (Int_t ibin=0;ibin<nptbins;ibin++){
+    for (Int_t ivar = 0; ivar<nvars; ivar++){
+      cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
+    }
+  }
+  
+  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+  SetUseSpecialCuts(kTRUE);
+  SetRemoveDaughtersFromPrim(kTRUE);
+  
+  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+  delete [] cutsMatrixTransposeStand;
+  cutsMatrixTransposeStand=NULL;
+
+  // PID SETTINGS
+  AliAODPidHF* pidObj=new AliAODPidHF();
+  //pidObj->SetName("pid4D0");
+  Int_t mode=1;
+  const Int_t nlims=2;
+  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+  Bool_t compat=kTRUE; //effective only for this mode
+  Bool_t asym=kTRUE;
+  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+  pidObj->SetMatch(mode);
+  pidObj->SetPLimit(plims,nlims);
+  pidObj->SetSigma(sigmas);
+  pidObj->SetCompat(compat);
+  pidObj->SetTPC(kTRUE);
+  pidObj->SetTOF(kTRUE);
+  
+  SetPidHF(pidObj);
+  SetUsePID(kTRUE);
+  SetUseDefaultPID(kFALSE);
+
+
+  PrintAll();
+
+  delete pidObj;
+  pidObj=NULL;
+
+  return;
+
+}
+
+
+void AliRDHFCutsD0toKpi::SetStandardCutsPbPb2010() {
+  //
+  //PRELIMINARY CUTS USED FOR 2010 PbPb analysis
+  //... EVOLVING SOON 
+  // 
+  
+  SetName("D0toKpiCutsStandard");
+  SetTitle("Standard Cuts for D0 analysis in PbPb2010 run");
+  
+  // PILE UP REJECTION
+  //SetOptPileup(AliRDHFCuts::kRejectPileupEvent);
+  // CENTRALITY SELECTION
+  SetMinCentrality(0.);
+  SetMaxCentrality(80.);
+  SetUseCentrality(AliRDHFCuts::kCentV0M);
+
+
+  // EVENT CUTS
+  SetMinVtxContr(1);
+  // MAX Z-VERTEX CUT
+  SetMaxVtxZ(10.);
+  
+  // TRACKS ON SINGLE TRACKS
+  AliESDtrackCuts *esdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts","default");
+  esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
+  esdTrackCuts->SetRequireTPCRefit(kTRUE);
+  esdTrackCuts->SetRequireITSRefit(kTRUE);
+  //  esdTrackCuts->SetMinNClustersITS(4);
+  esdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+  esdTrackCuts->SetMinDCAToVertexXY(0.);
+  esdTrackCuts->SetEtaRange(-0.8,0.8);
+  esdTrackCuts->SetPtRange(0.8,1.e10);
+
+  esdTrackCuts->SetMaxDCAToVertexXY(1.);  
+  esdTrackCuts->SetMaxDCAToVertexZ(1.);
+  esdTrackCuts->SetMinDCAToVertexXYPtDep("0.0100*TMath::Max(0.,(1-TMath::Floor(TMath::Abs(pt)/2.)))");  
+
+
+  AddTrackCuts(esdTrackCuts);
+
+  // CANDIDATE CUTS  
+  const Int_t nptbins =13;
+  const Double_t ptmax = 9999.;
+  const Int_t nvars=11;
+  Float_t ptbins[nptbins+1];
+  ptbins[0]=0.;
+  ptbins[1]=0.5;       
+  ptbins[2]=1.;
+  ptbins[3]=2.;
+  ptbins[4]=3.;
+  ptbins[5]=4.;
+  ptbins[6]=5.;
+  ptbins[7]=6.;
+  ptbins[8]=8.;
+  ptbins[9]=12.;
+  ptbins[10]=16.;
+  ptbins[11]=20.;
+  ptbins[12]=24.;
+  ptbins[13]=ptmax;
+
+  SetGlobalIndex(nvars,nptbins);
+  SetPtBins(nptbins+1,ptbins);
+  SetMinPtCandidate(2.);
+
+  Float_t cutsMatrixD0toKpiStand[nptbins][nvars]={{0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-50000.*1E-8,0.85,0.,5.},/* pt<0.5*/
+                                                 {0.400,400.*1E-4,0.8,0.3,0.3,1000.*1E-4,1000.*1E-4,-50000.*1E-8,0.85,0.,5.},/* 0.5<pt<1*/
+                                                 {0.400,400.*1E-4,0.8,0.4,0.4,1000.*1E-4,1000.*1E-4,-43000.*1E-8,0.85,0.,5.},/* 1<pt<2 */
+                                                 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-40000.*1E-8,0.95,0.998,5.},/* 2<pt<3 */
+                                                 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-36000.*1E-8,0.95,0.998,5.},/* 3<pt<4 */
+                                                 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-27000.*1E-8,0.95,0.998,5.},/* 4<pt<5 */
+                                                 {0.400,250.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-21000.*1E-8,0.92,0.998,5.},/* 5<pt<6 */
+                                                 {0.400,270.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-14000.*1E-8,0.88,0.998,5.},/* 6<pt<8 */
+                                                 {0.400,300.*1E-4,0.8,0.7,0.7,1000.*1E-4,1000.*1E-4,-5000.*1E-8,0.85,0.998,5.},/* 8<pt<12 */
+                                                 {0.400,350.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.83,0.998,5.},/* 12<pt<16 */
+                                                 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.82,0.998,5.},/* 16<pt<20 */
+                                                 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.81,0.998,5.},/* 20<pt<24 */
+                                                 {0.400,400.*1E-4,1.0,0.7,0.7,1000.*1E-4,1000.*1E-4,-1000.*1E-8,0.8,0.998,5.}};/* pt>24 */
+  
+  
+  //CREATE TRANSPOSE MATRIX...REVERSE INDICES as required by AliRDHFCuts
+  Float_t **cutsMatrixTransposeStand=new Float_t*[nvars];
+  for(Int_t iv=0;iv<nvars;iv++)cutsMatrixTransposeStand[iv]=new Float_t[nptbins];
+  
+  for (Int_t ibin=0;ibin<nptbins;ibin++){
+    for (Int_t ivar = 0; ivar<nvars; ivar++){
+      cutsMatrixTransposeStand[ivar][ibin]=cutsMatrixD0toKpiStand[ibin][ivar];      
+    }
+  }
+  
+  SetCuts(nvars,nptbins,cutsMatrixTransposeStand);
+  SetUseSpecialCuts(kTRUE);
+  SetRemoveDaughtersFromPrim(kFALSE);// THIS IS VERY IMPORTANT! TOO SLOW IN PbPb
+  for(Int_t iv=0;iv<nvars;iv++) delete [] cutsMatrixTransposeStand[iv];
+  delete [] cutsMatrixTransposeStand;
+  cutsMatrixTransposeStand=NULL;
+  
+  // PID SETTINGS
+  AliAODPidHF* pidObj=new AliAODPidHF();
+  //pidObj->SetName("pid4D0");
+  Int_t mode=1;
+  const Int_t nlims=2;
+  Double_t plims[nlims]={0.6,0.8}; //TPC limits in momentum [GeV/c]
+  Bool_t compat=kTRUE; //effective only for this mode
+  Bool_t asym=kTRUE;
+  Double_t sigmas[5]={2.,1.,0.,3.,0.}; //to be checked and to be modified with new implementation of setters by Rossella
+  pidObj->SetAsym(asym);// if you want to use the asymmetric bands in TPC
+  pidObj->SetMatch(mode);
+  pidObj->SetPLimit(plims,nlims);
+  pidObj->SetSigma(sigmas);
+  pidObj->SetCompat(compat);
+  pidObj->SetTPC(kTRUE);
+  pidObj->SetTOF(kTRUE);
+  
+  SetPidHF(pidObj);
+  SetUsePID(kTRUE);
+  SetUseDefaultPID(kFALSE);
+
+
+  PrintAll();
+
+
+  delete pidObj;
+  pidObj=NULL;
+
+  return;
+
+}