]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisVertexingHF.cxx
Updates D* cuts for peripheral PbPb events
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisVertexingHF.cxx
index 7f7a7500b99a44fbe45c4158938564f55c781f18..9c74c00b8ef5f680e0641681bd2487e25951de4a 100644 (file)
@@ -13,6 +13,8 @@
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
+/* $Id$ */
+
 //----------------------------------------------------------------------------
 //    Implementation of the heavy-flavour vertexing analysis class
 // Candidates are stored in the AOD as objects deriving from AliAODRecoDecay.
@@ -98,7 +100,9 @@ fFindVertexForCascades(kTRUE),
 fMassCutBeforeVertexing(kFALSE),
 fMassCalc2(0),
 fMassCalc3(0),
-fMassCalc4(0)
+fMassCalc4(0),
+fnTrksTotal(0),
+fnSeleTrksTotal(0)
 {
   // Default constructor
 
@@ -143,7 +147,9 @@ fFindVertexForCascades(source.fFindVertexForCascades),
 fMassCutBeforeVertexing(source.fMassCutBeforeVertexing),
 fMassCalc2(source.fMassCalc2),
 fMassCalc3(source.fMassCalc3),
-fMassCalc4(source.fMassCalc4)
+fMassCalc4(source.fMassCalc4),
+fnTrksTotal(0),
+fnSeleTrksTotal(0)
 {
   //
   // Copy constructor
@@ -205,7 +211,7 @@ AliAnalysisVertexingHF::~AliAnalysisVertexingHF() {
   if(fCutsLctoV0) { delete fCutsLctoV0; fCutsLctoV0=0; }
   if(fCutsD0toKpipipi) { delete fCutsD0toKpipipi; fCutsD0toKpipipi=0; }
   if(fCutsDStartoKpipi) { delete fCutsDStartoKpipi; fCutsDStartoKpipi=0; }
-  if(fAODMap) { delete fAODMap; fAODMap=0; }
+  if(fAODMap) { delete [] fAODMap; fAODMap=0; }
   if(fMassCalc2) { delete fMassCalc2; fMassCalc2=0; }
   if(fMassCalc3) { delete fMassCalc3; fMassCalc3=0; }
   if(fMassCalc4) { delete fMassCalc4; fMassCalc4=0; }
@@ -401,6 +407,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 
   trkEntries = (Int_t)event->GetNumberOfTracks();
   AliDebug(1,Form(" Number of tracks: %d",trkEntries));
+  fnTrksTotal += trkEntries;
 
   nv0 = (Int_t)event->GetNumberOfV0s();
   AliDebug(1,Form(" Number of V0s: %d",nv0));
@@ -422,8 +429,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   Int_t *evtNumber    = new Int_t[trkEntries];
   SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
     
-  //AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
-  printf(" Selected tracks: %d\n",nSeleTrks);
+  AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
+  fnSeleTrksTotal += nSeleTrks;
     
   TObjArray *twoTrackArray1    = new TObjArray(2);
   TObjArray *twoTrackArray2    = new TObjArray(2);
@@ -445,7 +452,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   // LOOP ON  POSITIVE  TRACKS
   for(iTrkP1=0; iTrkP1<nSeleTrks; iTrkP1++) {
 
-    if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));  
+    //if(iTrkP1%1==0) AliDebug(1,Form("  1st loop on pos: track number %d of %d",iTrkP1,nSeleTrks));  
     //if(iTrkP1%1==0) printf("  1st loop on pos: track number %d of %d\n",iTrkP1,nSeleTrks);  
 
     // get track from tracks array
@@ -459,7 +466,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
       // loop on V0's
       for(iv0=0; iv0<nv0; iv0++){
 
-       AliDebug(1,Form("   loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));     
+       //AliDebug(1,Form("   loop on v0s for track number %d and v0 number %d",iTrkP1,iv0));   
 
        // Get the V0 
        if(fInputAOD) {
@@ -580,7 +587,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        // Create and store the Cascade if passed the cuts
        ioCascade = MakeCascade(twoTrackArrayCasc,event,vertexCasc,v0,dcaCasc,okCascades);
        if(okCascades && ioCascade) {
-         AliDebug(1,Form("Storing a cascade object... "));
+         //AliDebug(1,Form("Storing a cascade object... "));
          // add the vertex and the cascade to the AOD
          AliAODVertex *vCasc = new(verticesHFRef[iVerticesHF++])AliAODVertex(*vertexCasc);
          rc = new(aodCascadesRef[iCascades++])AliAODRecoCascadeHF(*ioCascade);
@@ -616,7 +623,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
     // LOOP ON  NEGATIVE  TRACKS
     for(iTrkN1=0; iTrkN1<nSeleTrks; iTrkN1++) {
 
-      if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));  
+      //if(iTrkN1%1==0) AliDebug(1,Form("    1st loop on neg: track number %d of %d",iTrkN1,nSeleTrks));  
       //if(iTrkN1%1==0) printf("    1st loop on neg: track number %d of %d\n",iTrkN1,nSeleTrks);  
 
       if(iTrkN1==iTrkP1) continue;
@@ -690,6 +697,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            v2Prong->SetParent(rd);
            AddRefs(v2Prong,rd,event,twoTrackArray1);
          }
+         // Set selection bit for PID
+         if(okD0) SetSelectionBitForPID(fCutsD0toKpi,rd,AliRDHFCuts::kD0toKpiPID);
        }
        // D* candidates
        if(fDstar && okD0fromDstar && !isLikeSign2Prong) {
@@ -718,7 +727,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
                 evtNumber[iTrkP1]==evtNumber[iTrkN1]) continue;
            }
 
-           if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));  
+           //if(iTrkSoftPi%1==0) AliDebug(1,Form("    1st loop on pi_s: track number %d of %d",iTrkSoftPi,nSeleTrks));  
 
            trackD0->PropagateToDCA(fV1,fBzkG,kVeryBig);
            if(trackD0->GetSigmaY2()<0. || trackD0->GetSigmaZ2()<0.) continue; // this is insipired by the AliITStrackV2::Invariant() checks
@@ -771,6 +780,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
              if(!fInputAOD) vCasc->AddDaughter(rd); // just to fill ref #0 
              AddRefs(vCasc,rc,event,twoTrackArrayCasc);
              vCasc->AddDaughter(rd); // add the D0 (in ref #1)
+             // Set selection bit for PID
+             SetSelectionBitForPID(fCutsDStartoKpipi,rc,AliRDHFCuts::kDstarPID);
            }
            twoTrackArrayCasc->Clear();
            trackPi=0; 
@@ -798,7 +809,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 
        if(iTrkP2==iTrkP1 || iTrkP2==iTrkN1) continue;
 
-       if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));  
+       //if(iTrkP2%1==0) AliDebug(1,Form("    2nd loop on pos: track number %d of %d",iTrkP2,nSeleTrks));  
 
        // get track from tracks array
        postrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkP2);
@@ -852,7 +863,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            threeTrackArray->AddAt(postrack2,2);
          }
          if(fMassCutBeforeVertexing)
-           massCutOK = SelectInvMass(threeTrackArray);
+           massCutOK = SelectInvMassAndPt(threeTrackArray);
        }
 
        if(f3Prong && !massCutOK) {
@@ -891,6 +902,10 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
              v3Prong->SetParent(rd);
              AddRefs(v3Prong,rd,event,threeTrackArray);
            }
+           // Set selection bit for PID
+           SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
+           SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
+           SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
          }
          if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
          if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;} 
@@ -919,7 +934,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 
            if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 
-           if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
+           //if(iTrkN2%1==0) AliDebug(1,Form("    3rd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
 
            // get track from tracks array
            negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
@@ -955,7 +970,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            // check invariant mass cuts for D0
            massCutOK=kTRUE;
            if(fMassCutBeforeVertexing) 
-             massCutOK = SelectInvMass(fourTrackArray);
+             massCutOK = SelectInvMassAndPt(fourTrackArray);
            
            if(!massCutOK) {
              fourTrackArray->Clear(); 
@@ -998,7 +1013,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
 
        if(iTrkN2==iTrkP1 || iTrkN2==iTrkP2 || iTrkN2==iTrkN1) continue;
 
-       if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
+       //if(iTrkN2%1==0) AliDebug(1,Form("    2nd loop on neg: track number %d of %d",iTrkN2,nSeleTrks));  
 
        // get track from tracks array
        negtrack2 = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkN2);
@@ -1041,7 +1056,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        // check invariant mass cuts for D+,Ds,Lc
         massCutOK=kTRUE;
        if(fMassCutBeforeVertexing && f3Prong) 
-         massCutOK = SelectInvMass(threeTrackArray);
+         massCutOK = SelectInvMassAndPt(threeTrackArray);
 
        if(!massCutOK) { 
          threeTrackArray->Clear();
@@ -1076,6 +1091,10 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
              v3Prong->SetParent(rd);
              AddRefs(v3Prong,rd,event,threeTrackArray);
            }
+           // Set selection bit for PID
+           SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
+           SetSelectionBitForPID(fCutsDstoKKpi,rd,AliRDHFCuts::kDsPID);
+           SetSelectionBitForPID(fCutsLctopKpi,rd,AliRDHFCuts::kLcPID);
          }
          if(io3Prong) {delete io3Prong; io3Prong=NULL;} 
          if(secVert3PrAOD) {delete secVert3PrAOD; secVert3PrAOD=NULL;}
@@ -1147,6 +1166,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
     if(fAODMap) { delete fAODMap; fAODMap=NULL; }
   }
 
+  //printf("Trks: total %d  sele %d\n",fnTrksTotal,fnSeleTrksTotal);
+
   return;
 }
 //----------------------------------------------------------------------------
@@ -1203,6 +1224,78 @@ void AliAnalysisVertexingHF::AddDaughterRefs(AliAODVertex *v,
 
   return;
 }      
+//---------------------------------------------------------------------------
+void AliAnalysisVertexingHF::FixReferences(AliAODEvent *aod)  
+{
+  // Checks that the references to the daughter tracks are properly
+  // assigned and reassigns them if needed
+  //
+
+
+  TClonesArray *inputArray=(TClonesArray*)aod->GetList()->FindObject("VerticesHF");
+  if(!inputArray) return;
+
+  AliAODTrack *track = 0;
+  AliAODVertex *vertex = 0;
+
+  Bool_t needtofix=kFALSE;
+  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
+    vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
+    for(Int_t id=0; id<vertex->GetNDaughters(); id++) {
+      track = (AliAODTrack*)vertex->GetDaughter(id);
+      if(!track->GetStatus()) needtofix=kTRUE;
+    }
+    if(needtofix) break;
+  }
+
+  if(!needtofix) return;
+
+
+  printf("Fixing references\n");
+
+  fAODMapSize = 100000;
+  fAODMap = new Int_t[fAODMapSize];
+
+  for(Int_t i=0; i<aod->GetNumberOfTracks(); i++) {
+    track = aod->GetTrack(i);
+
+    // skip pure ITS SA tracks
+    if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
+
+    // skip tracks without ITS
+    if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
+
+    // TEMPORARY: check that the cov matrix is there
+    Double_t covtest[21];
+    if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
+    //
+
+    fAODMap[(Int_t)track->GetID()] = i;
+  }
+
+
+  Int_t ids[4]={-1,-1,-1,-1};
+  for(Int_t iv=0; iv<inputArray->GetEntriesFast(); iv++) {
+    Bool_t cascade=kFALSE;
+    vertex = (AliAODVertex*)inputArray->UncheckedAt(iv);
+    Int_t id=0;
+    Int_t nDgs = vertex->GetNDaughters();
+    for(id=0; id<nDgs; id++) {
+      track = (AliAODTrack*)vertex->GetDaughter(id);
+      if(track->Charge()==0) {cascade=kTRUE; continue;} // cascade
+      ids[id]=(Int_t)track->GetID();
+      vertex->RemoveDaughter(track);
+    }
+    if(cascade) continue;
+    for(id=0; id<nDgs; id++) {
+      track = aod->GetTrack(fAODMap[ids[id]]);
+      vertex->AddDaughter(track);
+    }
+    
+  }
+
+  return;
+}
 //----------------------------------------------------------------------------
 AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
                                   TObjArray *twoTrackArray,AliVEvent *event,
@@ -1244,6 +1337,7 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
   // select D*->D0pi
   if(fDstar) {
     okDstar = (Bool_t)fCutsDStartoKpipi->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
+    if(okDstar) theCascade->SetSelectionBit(AliRDHFCuts::kDstarCuts);
   }
   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
   tmpCascade->UnsetOwnPrimaryVtx(); 
@@ -1253,6 +1347,7 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
   }
   if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
   //---
+
   
   return theCascade;
 }
@@ -1299,14 +1394,13 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
   }
 
   // select Cascades
-  Bool_t okLcksp=0, okLcLpi=0;
   if(fCascades && fInputAOD){
     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
-    if(okCascades==1) okLcksp=1;
-    if(okCascades==2) okLcLpi=1;
-    if(okCascades==3) { okLcksp=1; okLcLpi=1;}
   }
-  else { AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); okCascades=kTRUE; }// no cuts implemented from ESDs
+  else { 
+    //AliDebug(2,Form("The cascade is contructed from ESDs, no cuts are applied")); 
+    okCascades=kTRUE; 
+  }// no cuts implemented from ESDs
   tmpCascade->GetSecondaryVtx()->RemoveDaughters();
   tmpCascade->UnsetOwnPrimaryVtx(); 
   delete tmpCascade; tmpCascade=NULL;
@@ -1346,11 +1440,11 @@ AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 
   // invariant mass cut (try to improve coding here..)
   Bool_t okMassCut=kFALSE;
-  if(!okMassCut && fD0toKpi)   if(SelectInvMass(0,2,px,py,pz)) okMassCut=kTRUE;
-  if(!okMassCut && fJPSItoEle) if(SelectInvMass(1,2,px,py,pz)) okMassCut=kTRUE;
-  if(!okMassCut && fDstar)     if(SelectInvMass(3,2,px,py,pz)) okMassCut=kTRUE;
+  if(!okMassCut && fD0toKpi)   if(SelectInvMassAndPt(0,2,px,py,pz)) okMassCut=kTRUE;
+  if(!okMassCut && fJPSItoEle) if(SelectInvMassAndPt(1,2,px,py,pz)) okMassCut=kTRUE;
+  if(!okMassCut && fDstar)     if(SelectInvMassAndPt(3,2,px,py,pz)) okMassCut=kTRUE;
   if(!okMassCut) {
-    AliDebug(2," candidate didn't pass mass cut");
+    //AliDebug(2," candidate didn't pass mass cut");
     return 0x0;    
   }
   // primary vertex to be used by this candidate
@@ -1375,14 +1469,26 @@ AliAODRecoDecayHF2Prong *AliAnalysisVertexingHF::Make2Prong(
 
  
   if(postrack->Charge()!=0 && negtrack->Charge()!=0) { // don't apply these cuts if it's a Dstar 
+
+    // Add daughter references already here
+    if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,twoTrackArray);
+
     // select D0->Kpi
-    if(fD0toKpi)   okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
+    if(fD0toKpi)   {
+      okD0 = (Bool_t)fCutsD0toKpi->IsSelected(the2Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event);
+      if(okD0) the2Prong->SetSelectionBit(AliRDHFCuts::kD0toKpiCuts);
+    }
     //if(fDebug && fD0toKpi) printf("   %d\n",(Int_t)okD0);
     // select J/psi from B
-    if(fJPSItoEle)   okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
+    if(fJPSItoEle)   {
+      okJPSI = (Bool_t)fCutsJpsitoee->IsSelected(the2Prong,AliRDHFCuts::kCandidate);
+    }
     //if(fDebug && fJPSItoEle) printf("   %d\n",(Int_t)okJPSI);
     // select D0->Kpi from Dstar
-    if(fDstar)   okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
+    if(fDstar)   {
+      okD0fromDstar = (Bool_t)fCutsDStartoKpipi->IsD0FromDStarSelected(the2Prong->Pt(),the2Prong,AliRDHFCuts::kCandidate);
+      if(okD0fromDstar) the2Prong->SetSelectionBit(AliRDHFCuts::kD0fromDstarCuts);
+    }
     //if(fDebug && fDstar) printf("   %d\n",(Int_t)okD0fromDstar);
   }
 
@@ -1443,9 +1549,9 @@ AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
   // invariant mass cut for D+, Ds, Lc
   Bool_t okMassCut=kFALSE;
   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed 
-  if(!okMassCut && f3Prong) if(SelectInvMass(2,3,px,py,pz)) okMassCut=kTRUE;
+  if(!okMassCut && f3Prong) if(SelectInvMassAndPt(2,3,px,py,pz)) okMassCut=kTRUE;
   if(!okMassCut) {
-    AliDebug(2," candidate didn't pass mass cut");
+    //AliDebug(2," candidate didn't pass mass cut");
     return 0x0;    
   }
 
@@ -1478,14 +1584,25 @@ AliAODRecoDecayHF3Prong* AliAnalysisVertexingHF::Make3Prong(
 
   delete primVertexAOD; primVertexAOD=NULL;
 
+  // Add daughter references already here
+  if(fInputAOD) AddDaughterRefs(secVert,(AliAODEvent*)event,threeTrackArray);
+
   // select D+->Kpipi, Ds->KKpi, Lc->pKpi
   if(f3Prong) {
     ok3Prong = kFALSE;
     
-    if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
-    if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
-    if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) ok3Prong = kTRUE;
-    
+    if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
+      ok3Prong = kTRUE;
+      the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
+    }
+    if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
+      ok3Prong = kTRUE;
+      the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
+    }
+    if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
+      ok3Prong = kTRUE;
+      the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
+    } 
   }
   //if(fDebug) printf("ok3Prong: %d\n",(Int_t)ok3Prong);
 
@@ -1554,7 +1671,7 @@ AliAODRecoDecayHF4Prong* AliAnalysisVertexingHF::Make4Prong(
   Bool_t okMassCut=kFALSE;
   if(fMassCutBeforeVertexing) okMassCut=kTRUE; // mass cut already done and passed 
   if(!okMassCut && !(fCutsD0toKpipipi->GetUsePID())) {      //no PID, to be implemented with PID
-    if(SelectInvMass(4,4,px,py,pz)) okMassCut=kTRUE;
+    if(SelectInvMassAndPt(4,4,px,py,pz)) okMassCut=kTRUE;
   }
   if(!okMassCut) {
     //if(fDebug) printf(" candidate didn't pass mass cut\n");
@@ -1711,7 +1828,7 @@ AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
 
     if(!vertexESD) return vertexAOD;
     if(vertexESD->GetNContributors()<=0) { 
-      AliDebug(2,"vertexing failed"); 
+      //AliDebug(2,"vertexing failed"); 
       delete vertexESD; vertexESD=NULL;
       return vertexAOD;
     }
@@ -1802,7 +1919,7 @@ AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkA
     if(!vertexESD) return vertexAOD;
 
     if(vertexESD->GetNContributors()!=trkArray->GetEntriesFast()) { 
-      AliDebug(2,"vertexing failed"); 
+      //AliDebug(2,"vertexing failed"); 
       delete vertexESD; vertexESD=NULL;
       return vertexAOD;
     }
@@ -1840,7 +1957,7 @@ AliAODVertex* AliAnalysisVertexingHF::ReconstructSecondaryVertex(TObjArray *trkA
   return vertexAOD;
 }
 //-----------------------------------------------------------------------------
-Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
+Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(TObjArray *trkArray) const {
   // Invariant mass cut on tracks
 
   Int_t nProngs=trkArray->GetEntriesFast();
@@ -1866,7 +1983,7 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
     px3[1] = momentum[0]; py3[1] = momentum[1]; pz3[1] = momentum[2]; 
     postrack2->GetPxPyPz(momentum);
     px3[2] = momentum[0]; py3[2] = momentum[1]; pz3[2] = momentum[2]; 
-    retval = SelectInvMass(2,3,px3,py3,pz3);
+    retval = SelectInvMassAndPt(2,3,px3,py3,pz3);
     break;
   case 4:
     // invariant mass cut for D0->4prong
@@ -1882,31 +1999,37 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(TObjArray *trkArray) const {
     px4[2] = momentum[0]; py4[2] = momentum[1]; pz4[2] = momentum[2];
     negtrack2->GetPxPyPz(momentum);
     px4[3] = momentum[0]; py4[3] = momentum[1]; pz4[3] = momentum[2];
-    retval = SelectInvMass(4,4,px4,py4,pz4);
+    retval = SelectInvMassAndPt(4,4,px4,py4,pz4);
     break;
   default:
-    printf("SelectInvMass(): wrong decay selection\n");
+    printf("SelectInvMassAndPt(): wrong decay selection\n");
     break;
   }
 
   return retval;
 }
 //-----------------------------------------------------------------------------
-Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
+Bool_t AliAnalysisVertexingHF::SelectInvMassAndPt(Int_t decay,
                                             Int_t nprongs,
                                             Double_t *px,
                                             Double_t *py,
                                             Double_t *pz) const {
-  // Check invariant mass cut
+  // Check invariant mass cut and pt candidate cut
 
   UInt_t pdg2[2],pdg3[3],pdg4[4];
   Double_t mPDG,minv;
+  Double_t minPt=0;
 
   Bool_t retval=kFALSE;
   switch (decay) 
     { 
     case 0:                  // D0->Kpi
       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
+      // pt cut
+      minPt=fCutsD0toKpi->GetMinPtCandidate();
+      if(minPt>0.1) 
+       if(fMassCalc2->Pt2() < minPt*minPt) break;
+      // mass cut
       pdg2[0]=211; pdg2[1]=321;
       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
       minv = fMassCalc2->InvMass(nprongs,pdg2);
@@ -1917,6 +2040,11 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       break;
     case 1:                  // JPSI->ee
       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
+      // pt cut
+      minPt=fCutsJpsitoee->GetMinPtCandidate();
+      if(minPt>0.1) 
+       if(fMassCalc2->Pt2() < minPt*minPt) break;
+      // mass cut
       pdg2[0]=11; pdg2[1]=11;
       mPDG=TDatabasePDG::Instance()->GetParticle(443)->Mass();
       minv = fMassCalc2->InvMass(nprongs,pdg2);
@@ -1924,6 +2052,12 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       break;
     case 2:                  // D+->Kpipi
       fMassCalc3->SetPxPyPzProngs(nprongs,px,py,pz);
+      // pt cut
+      minPt=TMath::Min(fCutsDplustoKpipi->GetMinPtCandidate(),fCutsDstoKKpi->GetMinPtCandidate());
+      minPt=TMath::Min(minPt,fCutsLctopKpi->GetMinPtCandidate());
+      if(minPt>0.1) 
+       if(fMassCalc3->Pt2() < minPt*minPt) break;
+      // mass cut
       pdg3[0]=211; pdg3[1]=321; pdg3[2]=211;
       mPDG=TDatabasePDG::Instance()->GetParticle(411)->Mass();
       minv = fMassCalc3->InvMass(nprongs,pdg3);
@@ -1947,6 +2081,11 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       break;
     case 3:                  // D*->D0pi
       fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
+      // pt cut
+      minPt=fCutsDStartoKpipi->GetMinPtCandidate();
+      if(minPt>0.1) 
+       if(fMassCalc2->Pt2() < minPt*minPt) break;
+      // mass cut
       pdg2[0]=211; pdg2[1]=421; // in twoTrackArrayCasc we put the pion first
       mPDG=TDatabasePDG::Instance()->GetParticle(413)->Mass();
       minv = fMassCalc2->InvMass(nprongs,pdg2);
@@ -1954,6 +2093,11 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       break;
     case 4:                 // D0->Kpipipi without PID
       fMassCalc4->SetPxPyPzProngs(nprongs,px,py,pz);
+      // pt cut
+      minPt=fCutsD0toKpipipi->GetMinPtCandidate();
+      if(minPt>0.1) 
+       if(fMassCalc4->Pt2() < minPt*minPt) break;
+      // mass cut
       pdg4[0]=321; pdg4[1]=211; pdg4[2]=211; pdg4[3]=211;
       mPDG=TDatabasePDG::Instance()->GetParticle(421)->Mass();
       minv = fMassCalc4->InvMass(nprongs,pdg4);
@@ -1972,7 +2116,7 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
       break;
     default:
-      printf("SelectInvMass(): wrong decay selection\n");
+      printf("SelectInvMassAndPt(): wrong decay selection\n");
       break;
     }
 
@@ -2023,6 +2167,9 @@ void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
     // skip pure ITS SA tracks
     if(track->GetStatus()&AliESDtrack::kITSpureSA) continue;
 
+    // skip tracks without ITS
+    if(!(track->GetStatus()&AliESDtrack::kITSin)) continue;
+
     // TEMPORARY: check that the cov matrix is there
     Double_t covtest[21];
     if(!track->GetCovarianceXYZPxPyPz(covtest)) continue;
@@ -2098,6 +2245,20 @@ void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
   return;
 }
 //-----------------------------------------------------------------------------
+void AliAnalysisVertexingHF::SetSelectionBitForPID(AliRDHFCuts *cuts,AliAODRecoDecayHF *rd,Int_t bit) {
+  //
+  // Set the selection bit for PID
+  //
+  if(cuts->GetPidHF()) {
+    Bool_t usepid=cuts->GetIsUsePID();
+    cuts->SetUsePID(kTRUE);
+    if(cuts->IsSelectedPID(rd))
+      rd->SetSelectionBit(bit);
+    cuts->SetUsePID(usepid);
+  }
+  return;
+}
+//-----------------------------------------------------------------------------
 Bool_t AliAnalysisVertexingHF::SingleTrkCuts(AliESDtrack *trk,
                                             Bool_t &okDisplaced,Bool_t &okSoftPi) const 
 {
@@ -2145,7 +2306,10 @@ AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArr
   esdV0->PxPyPz(pxpypz);
   Double_t cv[21]; for(int i=0; i<21; i++) cv[i]=0;
   AliNeutralTrackParam *trackesdV0 = new AliNeutralTrackParam(xyz,pxpypz,cv,0);
-  if(!trackesdV0) return 0;
+  if(!trackesdV0) {
+    delete vertexV0;
+    return 0;
+  }
   Double_t d0z0[2],covd0z0[3];
   AliAODVertex *primVertexAOD = PrimaryVertex(); 
   trackesdV0->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
@@ -2156,6 +2320,8 @@ AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArr
   AliExternalTrackParam *negV0track = (AliExternalTrackParam*)twoTrackArrayV0->UncheckedAt(1);
   if( !posV0track || !negV0track) {
     if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
+    delete vertexV0;
+    delete primVertexAOD;
     return 0;
   }
   posV0track->PropagateToDCA(primVertexAOD,fBzkG,kVeryBig,d0z0,covd0z0);
@@ -2174,8 +2340,8 @@ AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArr
   AliAODv0 *aodV0 = new AliAODv0(vertexV0,dcaV0Daughters,dcaV0ToPrimVertex,pmom,nmom,dcaV0DaughterToPrimVertex);
   aodV0->SetOnFlyStatus(esdV0->GetOnFlyStatus());
 
-  if(trackesdV0) {delete trackesdV0; trackesdV0=NULL;}
-  if(primVertexAOD) {delete primVertexAOD; primVertexAOD=NULL;}
+  delete trackesdV0;
+  delete primVertexAOD;
 
   return aodV0;
 }