]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/vertexingHF/AliAnalysisVertexingHF.cxx
introduce the option to analyze with PROOF, plus cosmetics
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliAnalysisVertexingHF.cxx
index 061e105f1fef3b0a121be74a7011709d8370f159..bf3ef098394fe8b42aa47a88d0b1406e578eecc5 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.
@@ -81,6 +83,7 @@ f4Prong(kTRUE),
 fDstar(kTRUE),
 fCascades(kTRUE),
 fLikeSign(kFALSE),
+fLikeSign3prong(kFALSE),
 fMixEvent(kFALSE),
 fTrackFilter(0x0),
 fTrackFilterSoftPi(0x0),
@@ -98,11 +101,15 @@ fFindVertexForCascades(kTRUE),
 fMassCutBeforeVertexing(kFALSE),
 fMassCalc2(0),
 fMassCalc3(0),
-fMassCalc4(0)
+fMassCalc4(0),
+fnTrksTotal(0),
+fnSeleTrksTotal(0)
 {
   // Default constructor
 
-  Double_t d02[2],d03[3],d04[4];
+  Double_t d02[2]={0.,0.};
+  Double_t d03[3]={0.,0.,0.};
+  Double_t d04[4]={0.,0.,0.,0.};
   fMassCalc2 = new AliAODRecoDecay(0x0,2,0,d02);
   fMassCalc3 = new AliAODRecoDecay(0x0,3,1,d03);
   fMassCalc4 = new AliAODRecoDecay(0x0,4,0,d04);
@@ -126,6 +133,7 @@ f4Prong(source.f4Prong),
 fDstar(source.fDstar),
 fCascades(source.fCascades),
 fLikeSign(source.fLikeSign),
+fLikeSign3prong(source.fLikeSign3prong),
 fMixEvent(source.fMixEvent),
 fTrackFilter(source.fTrackFilter),
 fTrackFilterSoftPi(source.fTrackFilterSoftPi),
@@ -143,7 +151,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
@@ -170,6 +180,7 @@ AliAnalysisVertexingHF &AliAnalysisVertexingHF::operator=(const AliAnalysisVerte
   fDstar = source.fDstar;
   fCascades = source.fCascades;
   fLikeSign = source.fLikeSign;
+  fLikeSign3prong = source.fLikeSign3prong;
   fMixEvent = source.fMixEvent;
   fTrackFilter = source.fTrackFilter;
   fTrackFilterSoftPi = source.fTrackFilterSoftPi;
@@ -315,8 +326,8 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
     printf("ERROR: no aodLikeSign2ProngTClArr");
     return;
   }
-  if(fLikeSign && f3Prong && !aodLikeSign3ProngTClArr) {
-    printf("ERROR: no aodLikeSign2ProngTClArr");
+  if(fLikeSign3prong && f3Prong && !aodLikeSign3ProngTClArr) {
+    printf("ERROR: no aodLikeSign3ProngTClArr");
     return;
   }
 
@@ -353,7 +364,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
     aodLikeSign2ProngTClArr->Delete();                     
     iLikeSign2Prong = aodLikeSign2ProngTClArr->GetEntriesFast(); 
   }  
-  if(fLikeSign && f3Prong) {                                
+  if(fLikeSign3prong && f3Prong) {                                
     aodLikeSign3ProngTClArr->Delete();                     
     iLikeSign3Prong = aodLikeSign3ProngTClArr->GetEntriesFast(); 
   }  
@@ -401,6 +412,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));
@@ -417,20 +429,23 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   // for displaced tracks and soft pions (both charges) for D*,
   // and retrieves primary vertex
   TObjArray seleTrksArray(trkEntries);
+  TObjArray tracksAtVertex(trkEntries);
   UChar_t  *seleFlags = new UChar_t[trkEntries]; // bit 0: displaced, bit 1: softpi
   Int_t     nSeleTrks=0;
   Int_t *evtNumber    = new Int_t[trkEntries];
-  SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,nSeleTrks,seleFlags,evtNumber);
-    
+  SelectTracksAndCopyVertex(event,trkEntries,seleTrksArray,tracksAtVertex,nSeleTrks,seleFlags,evtNumber);
+
   AliDebug(1,Form(" Selected tracks: %d",nSeleTrks));
+  fnSeleTrksTotal += nSeleTrks;
     
+
   TObjArray *twoTrackArray1    = new TObjArray(2);
   TObjArray *twoTrackArray2    = new TObjArray(2);
   TObjArray *twoTrackArrayV0   = new TObjArray(2);
   TObjArray *twoTrackArrayCasc = new TObjArray(2);
   TObjArray *threeTrackArray   = new TObjArray(3);
   TObjArray *fourTrackArray    = new TObjArray(4);
-  
+
   Double_t dispersion;
   Bool_t isLikeSign2Prong=kFALSE,isLikeSign3Prong=kFALSE;
 
@@ -444,7 +459,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
@@ -458,7 +473,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) {
@@ -550,7 +565,6 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        // Fill in the object array to create the cascade
        twoTrackArrayCasc->AddAt(postrack1,0);
        twoTrackArrayCasc->AddAt(trackV0,1);
-
        // Compute the cascade vertex
        AliAODVertex *vertexCasc = 0;
        if(fFindVertexForCascades) {  
@@ -572,6 +586,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
          delete negV0track; negV0track=NULL;
          delete trackV0; trackV0=NULL;
          if(!fInputAOD) {delete v0; v0=NULL;}
+         twoTrackArrayV0->Clear();
          twoTrackArrayCasc->Clear();
          continue; 
        }
@@ -579,7 +594,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);
@@ -615,7 +630,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;
@@ -645,8 +660,10 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
       }
 
       // back to primary vertex
-      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+      //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+      //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+      SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
+      SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
 
       // DCA between the two tracks
       dcap1n1 = postrack1->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
@@ -656,7 +673,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
       twoTrackArray1->AddAt(postrack1,0);
       twoTrackArray1->AddAt(negtrack1,1);
       AliAODVertex *vertexp1n1 = ReconstructSecondaryVertex(twoTrackArray1,dispersion);
-      if(!vertexp1n1) { 
+      if(!vertexp1n1) {
        twoTrackArray1->Clear();
        negtrack1=0; 
        continue; 
@@ -719,15 +736,15 @@ 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
 
            // get track from tracks array
            trackPi = (AliESDtrack*)seleTrksArray.UncheckedAt(iTrkSoftPi);
-           trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
-
+           //      trackPi->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           SetParametersAtVertex(trackPi,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkSoftPi));
            twoTrackArrayCasc->AddAt(trackPi,0);
            twoTrackArrayCasc->AddAt(trackD0,1);
 
@@ -801,7 +818,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);
@@ -817,6 +834,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        }
 
        if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
+         if(!fLikeSign3prong) continue;
          if(postrack1->Charge()>0) { // ok: like-sign triplet (+++)
            isLikeSign3Prong=kTRUE;
          } else { // not ok
@@ -832,9 +850,13 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        }
 
        // back to primary vertex
-       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-       postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
-       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       //      postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
+       SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
+       SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
+      
        //printf("********** %d %d %d\n",postrack1->GetID(),postrack2->GetID(),negtrack1->GetID());
 
        dcap2n1 = postrack2->GetDCA(negtrack1,fBzkG,xdummy,ydummy);
@@ -855,7 +877,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            threeTrackArray->AddAt(postrack2,2);
          }
          if(fMassCutBeforeVertexing)
-           massCutOK = SelectInvMass(threeTrackArray);
+           massCutOK = SelectInvMassAndPt(threeTrackArray);
        }
 
        if(f3Prong && !massCutOK) {
@@ -889,10 +911,12 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
              v3Prong->SetParent(rd);
              AddRefs(v3Prong,rd,event,threeTrackArray);
            } else { // isLikeSign3Prong
-             rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
-             rd->SetSecondaryVtx(v3Prong);
-             v3Prong->SetParent(rd);
-             AddRefs(v3Prong,rd,event,threeTrackArray);
+             if(fLikeSign3prong){
+               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
+               rd->SetSecondaryVtx(v3Prong);
+               v3Prong->SetParent(rd);
+               AddRefs(v3Prong,rd,event,threeTrackArray);
+             }
            }
            // Set selection bit for PID
            SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
@@ -912,13 +936,17 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
           && dcap2n1 < fCutsD0toKpipipi->GetDCACut()) {
 
          // back to primary vertex
-         postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-         postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
-         negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         //      postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         //      postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         //      negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+         SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
+         SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
+         SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
+
          // Vertexing for these 3 (can be taken from above?)
           threeTrackArray->AddAt(postrack1,0);
           threeTrackArray->AddAt(negtrack1,1);
-          threeTrackArray->AddAt(postrack2,2);
+         threeTrackArray->AddAt(postrack2,2);
           AliAODVertex* vertexp1n1p2 = ReconstructSecondaryVertex(threeTrackArray,dispersion);
 
          // 3rd LOOP  ON  NEGATIVE  TRACKS (for 4 prong) 
@@ -926,7 +954,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);
@@ -944,10 +972,15 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
            }
 
            // back to primary vertex
-           postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-           postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
-           negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-           negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           // postrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+           SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
+           SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
+           SetParametersAtVertex(postrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP2));
+           SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
+
            dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
            if(dcap1n2 > fCutsD0toKpipipi->GetDCACut()) { negtrack2=0; continue; }
             dcap2n2 = postrack2->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
@@ -962,7 +995,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(); 
@@ -1005,7 +1038,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);
@@ -1021,6 +1054,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        }
 
        if(isLikeSign2Prong) { // like-sign pair -> have to build only like-sign triplet 
+         if(!fLikeSign3prong) continue;
          if(postrack1->Charge()<0) { // ok: like-sign triplet (---)
            isLikeSign3Prong=kTRUE;
          } else { // not ok
@@ -1031,9 +1065,12 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
        }
 
        // back to primary vertex
-       postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-       negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
-       negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       // postrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       // negtrack1->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       // negtrack2->PropagateToDCA(fV1,fBzkG,kVeryBig);
+       SetParametersAtVertex(postrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkP1));
+       SetParametersAtVertex(negtrack1,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN1));
+       SetParametersAtVertex(negtrack2,(AliExternalTrackParam*)tracksAtVertex.UncheckedAt(iTrkN2));
        //printf("********** %d %d %d\n",postrack1->GetID(),negtrack1->GetID(),negtrack2->GetID());
 
        dcap1n2 = postrack1->GetDCA(negtrack2,fBzkG,xdummy,ydummy);
@@ -1048,7 +1085,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();
@@ -1078,10 +1115,12 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
              v3Prong->SetParent(rd);
              AddRefs(v3Prong,rd,event,threeTrackArray);
            } else { // isLikeSign3Prong
-             rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
-             rd->SetSecondaryVtx(v3Prong);
-             v3Prong->SetParent(rd);
-             AddRefs(v3Prong,rd,event,threeTrackArray);
+             if(fLikeSign3prong){
+               rd = new(aodLikeSign3ProngRef[iLikeSign3Prong++])AliAODRecoDecayHF3Prong(*io3Prong);
+               rd->SetSecondaryVtx(v3Prong);
+               v3Prong->SetParent(rd);
+               AddRefs(v3Prong,rd,event,threeTrackArray);
+             }
            }
            // Set selection bit for PID
            SetSelectionBitForPID(fCutsDplustoKpipi,rd,AliRDHFCuts::kDplusPID);
@@ -1137,7 +1176,7 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
     AliDebug(1,Form(" Like-sign 2Prong in event = %d;\n",
                    (Int_t)aodLikeSign2ProngTClArr->GetEntriesFast()));
   }
-  if(fLikeSign && f3Prong) {
+  if(fLikeSign3prong && f3Prong) {
     AliDebug(1,Form(" Like-sign 3Prong in event = %d;\n",
                    (Int_t)aodLikeSign3ProngTClArr->GetEntriesFast()));
   }
@@ -1152,11 +1191,15 @@ void AliAnalysisVertexingHF::FindCandidates(AliVEvent *event,
   fourTrackArray->Delete();  delete fourTrackArray;
   delete [] seleFlags; seleFlags=NULL;
   if(evtNumber) {delete [] evtNumber; evtNumber=NULL;}
+  tracksAtVertex.Delete();
 
   if(fInputAOD) {
     seleTrksArray.Delete(); 
     if(fAODMap) { delete fAODMap; fAODMap=NULL; }
   }
+  
+
+  //printf("Trks: total %d  sele %d\n",fnTrksTotal,fnSeleTrksTotal);
 
   return;
 }
@@ -1315,8 +1358,15 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
   //--- selection cuts
   //
   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);
-  tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
+  if(fInputAOD){
+    Int_t idSoftPi=(Int_t)trackPi->GetID();
+    AliAODTrack* trackPiAOD=(AliAODTrack*)event->GetTrack(fAODMap[idSoftPi]);
+    tmpCascade->GetSecondaryVtx()->AddDaughter(trackPiAOD);
+  }else{
+    tmpCascade->GetSecondaryVtx()->AddDaughter(trackPi);
+  }
   tmpCascade->GetSecondaryVtx()->AddDaughter(rd2Prong);
+
   AliAODVertex *primVertexAOD=0;
   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
     // take event primary vertex
@@ -1372,9 +1422,17 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
 
   //--- selection cuts
   //
+
   AliAODRecoCascadeHF *tmpCascade = new AliAODRecoCascadeHF(*theCascade);  
-  tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
+  if(fInputAOD){
+    Int_t idBachelor=(Int_t)trackBachelor->GetID();
+    AliAODTrack* trackBachelorAOD=(AliAODTrack*)event->GetTrack(fAODMap[idBachelor]);
+    tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelorAOD);
+  }else{
+    tmpCascade->GetSecondaryVtx()->AddDaughter(trackBachelor);
+  }
   tmpCascade->GetSecondaryVtx()->AddDaughter(v0);
+  
   AliAODVertex *primVertexAOD=0;
   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
     // take event primary vertex
@@ -1382,12 +1440,15 @@ AliAODRecoCascadeHF* AliAnalysisVertexingHF::MakeCascade(
     if(!primVertexAOD) primVertexAOD = (AliAODVertex*)event->GetPrimaryVertex(); 
     tmpCascade->SetOwnPrimaryVtx(primVertexAOD);
   }
-
+  
   // select Cascades
   if(fCascades && fInputAOD){
     okCascades = (Bool_t)fCutsLctoV0->IsSelected(tmpCascade,AliRDHFCuts::kCandidate);
   }
-  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;
@@ -1427,11 +1488,12 @@ 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 && fCascades)  if(SelectInvMassAndPt(5,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
@@ -1456,9 +1518,13 @@ 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);
+      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);
@@ -1532,9 +1598,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;    
   }
 
@@ -1567,19 +1633,22 @@ 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)) {
+    if(fCutsDplustoKpipi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
       ok3Prong = kTRUE;
       the3Prong->SetSelectionBit(AliRDHFCuts::kDplusCuts);
     }
-    if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) {
+    if(fCutsDstoKKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
       ok3Prong = kTRUE;
       the3Prong->SetSelectionBit(AliRDHFCuts::kDsCuts);
     }
-    if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate)) {
+    if(fCutsLctopKpi->IsSelected(the3Prong,AliRDHFCuts::kCandidate,(AliAODEvent*)event)) {
       ok3Prong = kTRUE;
       the3Prong->SetSelectionBit(AliRDHFCuts::kLcCuts);
     } 
@@ -1651,7 +1720,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");
@@ -1778,11 +1847,12 @@ AliAODVertex* AliAnalysisVertexingHF::PrimaryVertex(const TObjArray *trkArray,
          }
        }
       }
+      for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
       //
       vertexer->SetSkipTracks(nTrksToSkip,skipped);
       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
       
-    } else if(fRmTrksFromPrimVtx) { 
+    } else if(fRmTrksFromPrimVtx && nTrks>0) { 
       // removing the prongs tracks
       
       TObjArray rmArray(nTrks);
@@ -1808,7 +1878,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;
     }
@@ -1899,7 +1969,14 @@ 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;
+    }
+    
+    Double_t vertRadius2=vertexESD->GetXv()*vertexESD->GetXv()+vertexESD->GetYv()*vertexESD->GetYv();
+    if(vertRadius2>8.){
+      // vertex outside beam pipe, reject candidate to avoid propagation through material
       delete vertexESD; vertexESD=NULL;
       return vertexAOD;
     }
@@ -1937,7 +2014,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();
@@ -1963,7 +2040,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
@@ -1979,31 +2056,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);
@@ -2014,6 +2097,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);
@@ -2021,6 +2109,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);
@@ -2044,6 +2138,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);
@@ -2051,6 +2150,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);
@@ -2068,8 +2172,18 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
       minv = fMassCalc4->InvMass(nprongs,pdg4);
       if(TMath::Abs(minv-mPDG)<fCutsD0toKpipipi->GetMassCut()) retval=kTRUE;
       break;
+    case 5:
+      fMassCalc2->SetPxPyPzProngs(nprongs,px,py,pz);
+      mPDG=TDatabasePDG::Instance()->GetParticle(4122)->Mass();
+      pdg2[0]=2212;pdg2[1]=310;
+      minv=fMassCalc2->InvMass(2,pdg2);
+      if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
+      pdg2[0]=211;pdg2[1]=3122;
+      minv=fMassCalc2->InvMass(2,pdg2);
+      if(TMath::Abs(minv-mPDG)<fCutsLctoV0->GetMassCut()) retval=kTRUE;
+      break;
     default:
-      printf("SelectInvMass(): wrong decay selection\n");
+      printf("SelectInvMassAndPt(): wrong decay selection\n");
       break;
     }
 
@@ -2078,8 +2192,9 @@ Bool_t AliAnalysisVertexingHF::SelectInvMass(Int_t decay,
 //-----------------------------------------------------------------------------
 void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
                                                       Int_t trkEntries,
-                                  TObjArray &seleTrksArray,Int_t &nSeleTrks,
-                                  UChar_t *seleFlags,Int_t *evtNumber)
+                                                      TObjArray &seleTrksArray,TObjArray &tracksAtVertex,
+                                                      Int_t &nSeleTrks,
+                                                      UChar_t *seleFlags,Int_t *evtNumber)
 {
   // Apply single-track preselection.
   // Fill a TObjArray with selected tracks (for displaced vertices or
@@ -2103,6 +2218,7 @@ void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
     vprimary->GetCovarianceMatrix(cov);
     fV1 = new AliESDVertex(pos,cov,100.,100,vprimary->GetName());
     indices = new UShort_t[event->GetNumberOfTracks()];
+    for(Int_t ijk=0; ijk<event->GetNumberOfTracks(); ijk++) indices[ijk]=0;
     fAODMapSize = 100000;
     fAODMap = new Int_t[fAODMapSize];
   }
@@ -2166,7 +2282,10 @@ void AliAnalysisVertexingHF::SelectTracksAndCopyVertex(const AliVEvent *event,
     }
 
     if(SingleTrkCuts(esdt,okDisplaced,okSoftPi) && nSeleTrks<trkEntries) {
+      esdt->PropagateToDCA(fV1,fBzkG,kVeryBig);
       seleTrksArray.AddLast(esdt);
+      tracksAtVertex.AddLast(new AliExternalTrackParam(*esdt));
+
       seleFlags[nSeleTrks]=0;
       if(okDisplaced) SETBIT(seleFlags[nSeleTrks],kBitDispl);
       if(okSoftPi)    SETBIT(seleFlags[nSeleTrks],kBitSoftPi);
@@ -2299,3 +2418,13 @@ AliAODv0* AliAnalysisVertexingHF::TransformESDv0toAODv0(AliESDv0 *esdV0, TObjArr
   return aodV0;
 }
 //-----------------------------------------------------------------------------
+void AliAnalysisVertexingHF::SetParametersAtVertex(AliESDtrack* esdt, AliExternalTrackParam* extpar) const{
+  // Set the stored track parameters at primary vertex into AliESDtrack
+  Double_t xyz[3], pxpypz[3], cv[21]; 
+  extpar->PxPyPz(pxpypz);                        
+  extpar->XvYvZv(xyz);
+  extpar->GetCovarianceXYZPxPyPz(cv);
+  Short_t sign=extpar->Charge();
+  esdt->Set(xyz,pxpypz,cv,sign);
+  return;
+}