]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEextraCuts.cxx
Update of hfe code
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEextraCuts.cxx
index f812754dcab404d72a8b6b65d3a930aa401d0c92..df8a6545060b5c24cfec9262c3921bc38b4b879c 100644 (file)
 #include <TString.h>
 #include <TMath.h>
 
+#include "AliAODEvent.h"
 #include "AliAODTrack.h"
 #include "AliAODPid.h"
+#include "AliAODVertex.h"
+#include "AliAODTrack.h"
 #include "AliESDEvent.h"
 #include "AliESDVertex.h"
 #include "AliESDtrack.h"
 #include "AliVParticle.h"
 #include "AliVertexerTracks.h"
 #include "AliVVertex.h"
+#include "AliExternalTrackParam.h"
 
 #include "AliHFEextraCuts.h"
 
 ClassImp(AliHFEextraCuts)
 
-const Int_t AliHFEextraCuts::fgkNQAhistos = 8;
+const Int_t AliHFEextraCuts::fgkNQAhistos = 10;
 
 //______________________________________________________
 AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
@@ -59,11 +63,17 @@ AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
   fMinNClustersTPCPID(0),
   fClusterRatioTPC(0.),
   fMinTrackletsTRD(0),
+  fMaxChi2TRD(5.0),
+  fMinNbITScls(0),
   fTRDtrackletsExact(0),
   fPixelITS(0),
+  fDriftITS(0),
   fTPCclusterDef(0),
   fTPCclusterRatioDef(0),
   fFractionTPCShared(-1.0),
+  fAbsHFEImpactParamNsigmaR(kTRUE),
+  fTOFsignalDx(1.0),
+  fTOFsignalDz(1.0),
   fCheck(kFALSE),
   fQAlist(0x0) ,
   fDebugLevel(0)
@@ -71,6 +81,7 @@ AliHFEextraCuts::AliHFEextraCuts(const Char_t *name, const Char_t *title):
   //
   // Default Constructor
   //
+  //printf("Set the number of min ITS clusters %d\n",(Int_t)fMinNbITScls);
   memset(fImpactParamCut, 0, sizeof(Float_t) * 4);
   memset(fIPcutParam, 0, sizeof(Float_t) * 4);
 }
@@ -85,11 +96,17 @@ AliHFEextraCuts::AliHFEextraCuts(const AliHFEextraCuts &c):
   fMinNClustersTPCPID(c.fMinNClustersTPCPID),
   fClusterRatioTPC(c.fClusterRatioTPC),
   fMinTrackletsTRD(c.fMinTrackletsTRD),
+  fMaxChi2TRD(c.fMaxChi2TRD),
+  fMinNbITScls(c.fMinNbITScls),
   fTRDtrackletsExact(c.fTRDtrackletsExact),
   fPixelITS(c.fPixelITS),
+  fDriftITS(c.fDriftITS),
   fTPCclusterDef(c.fTPCclusterDef),
   fTPCclusterRatioDef(c.fTPCclusterRatioDef),
   fFractionTPCShared(c.fFractionTPCShared),
+  fAbsHFEImpactParamNsigmaR(c.fAbsHFEImpactParamNsigmaR),
+  fTOFsignalDx(c.fTOFsignalDx),
+  fTOFsignalDz(c.fTOFsignalDz),
   fCheck(c.fCheck),
   fQAlist(0x0),
   fDebugLevel(0)
@@ -121,11 +138,17 @@ AliHFEextraCuts &AliHFEextraCuts::operator=(const AliHFEextraCuts &c){
     fMinNClustersTPC = c.fMinNClustersTPC;
     fMinNClustersTPCPID = c.fMinNClustersTPCPID;
     fMinTrackletsTRD = c.fMinTrackletsTRD;
+    fMaxChi2TRD      = c.fMaxChi2TRD;
+    fMinNbITScls = c.fMinNbITScls;
     fTRDtrackletsExact = c.fTRDtrackletsExact;
     fPixelITS = c.fPixelITS;
+    fDriftITS = c.fDriftITS;
     fTPCclusterDef = c.fTPCclusterDef;
     fTPCclusterRatioDef = c.fTPCclusterRatioDef;
     fFractionTPCShared = c.fFractionTPCShared;
+    fAbsHFEImpactParamNsigmaR = c.fAbsHFEImpactParamNsigmaR;
+    fTOFsignalDx = c.fTOFsignalDx;
+    fTOFsignalDz = c.fTOFsignalDz;
     fCheck = c.fCheck;
     fDebugLevel = c.fDebugLevel;
     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
@@ -170,9 +193,9 @@ Bool_t AliHFEextraCuts::IsSelected(TObject *o){
   //
   // Steering function for the track selection
   //
-  TString type = o->IsA()->GetName();
-  AliDebug(2, Form("Object type %s", type.Data()));
-  if(!type.CompareTo("AliESDtrack") || !type.CompareTo("AliAODTrack")){
+  TClass *type = o->IsA();
+  AliDebug(2, Form("Object type %s", type->GetName()));
+  if((type == AliESDtrack::Class()) || (type == AliAODTrack::Class())){
     return CheckRecCuts(dynamic_cast<AliVTrack *>(o));
   }
   return CheckMCCuts(dynamic_cast<AliVParticle *>(o));
@@ -186,11 +209,12 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
   // QA histograms are filled before track selection and for
   // selected tracks after track selection
   //
-  AliDebug(1, "Called");
+  AliDebug(1, Form("%s: Called", GetName()));
   ULong64_t survivedCut = 0;   // Bitmap for cuts which are passed by the track, later to be compared with fRequirements
   if(IsQAOn()) FillQAhistosRec(track, kBeforeCuts);
   // Apply cuts
-  Float_t impactR, impactZ;
+  Float_t impactR = -999.;
+  Float_t impactZ = -999.;
   Double_t hfeimpactR, hfeimpactnsigmaR;
   Double_t hfeimpactRcut, hfeimpactnsigmaRcut;
   Double_t maximpactRcut; 
@@ -200,18 +224,25 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
   }
+  Int_t  nclsITS = GetITSNbOfcls(track);
   UInt_t nclsTPC = GetTPCncls(track);
-  UInt_t nclsTPCPID = GetTPCnclusdEdx(track);
+  UInt_t nclsTPCPID = track->GetTPCsignalN();
   // printf("Check TPC findable clusters: %d, found Clusters: %d\n", track->GetTPCNclsF(), track->GetTPCNcls());
   Float_t fractionSharedClustersTPC = GetTPCsharedClustersRatio(track);
   Double_t ratioTPC = GetTPCclusterRatio(track);
   UChar_t trdTracklets;
   trdTracklets = track->GetTRDntrackletsPID();
+  Float_t trdchi2=-999.;
+  trdchi2=GetTRDchi(track);
   UChar_t itsPixel = track->GetITSClusterMap();
   Int_t status1 = GetITSstatus(track, 0);
   Int_t status2 = GetITSstatus(track, 1);
   Bool_t statusL0 = CheckITSstatus(status1);
   Bool_t statusL1 = CheckITSstatus(status2);
+  Double_t tofsignalDx = 0.0;
+  Double_t tofsignalDz = 0.0;
+  GetTOFsignalDxDz(track,tofsignalDx,tofsignalDz);
+
   if(TESTBIT(fRequirements, kTPCfractionShared)) {
     // cut on max fraction of shared TPC clusters
     if(TMath::Abs(fractionSharedClustersTPC) < fFractionTPCShared) SETBIT(survivedCut, kTPCfractionShared);    
@@ -243,7 +274,27 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
   }
   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
     // cut on max. HFE Impact Parameter n sigma in Radial direction
-    if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+    if(fAbsHFEImpactParamNsigmaR) {
+//      if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
+      if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
+        SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+//        printf("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+      }
+    }
+    else {
+      if(hfeimpactnsigmaRcut>0){
+        if(hfeimpactnsigmaR >= hfeimpactnsigmaRcut) {
+          SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+          //printf("1: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+        }
+      }
+      else{
+        if(hfeimpactnsigmaR <= hfeimpactnsigmaRcut) {
+          SETBIT(survivedCut, kMinHFEImpactParamNsigmaR);
+          //printf("2: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+        }
+      }
+    }
   }
   if(TESTBIT(fRequirements, kClusterRatioTPC)){
     // cut on min ratio of found TPC clusters vs findable TPC clusters
@@ -251,36 +302,66 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
   }
   if(TESTBIT(fRequirements, kMinTrackletsTRD)){
     // cut on minimum number of TRD tracklets
-    AliDebug(1, Form("Min TRD cut: [%d|%d], Require exact number [%s]\n", fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
+    AliDebug(1, Form("%s: Min TRD cut: [%d|%d], Require exact number [%s]\n", GetName(), fMinTrackletsTRD, trdTracklets, fTRDtrackletsExact ? "Yes" : "No"));
     if(fTRDtrackletsExact){
       if(trdTracklets == fMinTrackletsTRD) {
         SETBIT(survivedCut, kMinTrackletsTRD);
-        AliDebug(1, "Track Selected");
+        AliDebug(1, Form("%s: Track Selected", GetName()));
       }
     }else{
       if(trdTracklets >= fMinTrackletsTRD){ 
         SETBIT(survivedCut, kMinTrackletsTRD);
-        AliDebug(1, "Track Selected");
+        AliDebug(1, Form("%s: Track Selected", GetName()));
       }
       //printf("Min number of tracklets %d\n",fMinTrackletsTRD);
     }
   }
+
+  if(TESTBIT(fRequirements, kMaxTRDChi2)){
+    // cut on TRD chi2
+    AliDebug(1, Form("%s: Cut on TRD chi2: [%f|%f]\n", GetName(),fMaxChi2TRD, trdchi2));
+    if(trdchi2 < fMaxChi2TRD) {
+       SETBIT(survivedCut, kMaxTRDChi2);
+        AliDebug(1,Form("%s: survived %f\n",GetName(),trdchi2));
+    }
+  }
+
+  if(TESTBIT(fRequirements, kMinNbITScls)){
+    // cut on minimum number of ITS clusters
+    //printf(Form("Min ITS clusters: [%d|%d]\n", (Int_t)fMinNbITScls, nclsITS));
+    AliDebug(1, Form("%s: Min ITS clusters: [%d|%d]\n", GetName(), fMinNbITScls, nclsITS));
+    if(nclsITS >= ((Int_t)fMinNbITScls)) SETBIT(survivedCut, kMinNbITScls);
+  }
+  
   if(TESTBIT(fRequirements, kMinNClustersTPC)){
-    // cut on minimum number of TRD tracklets
-    AliDebug(1, Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
+    // cut on minimum number of TPC tracklets
+    //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
+    AliDebug(1, Form("%s: Min TPC cut: [%d|%d]\n", GetName(), fMinNClustersTPC, nclsTPC));
     if(nclsTPC >= fMinNClustersTPC) SETBIT(survivedCut, kMinNClustersTPC);
   }
   if(TESTBIT(fRequirements, kMinNClustersTPCPID)){
-    AliDebug(1, Form("Min TPC PID cut: [%d|%d]\n", fMinNClustersTPCPID, nclsTPCPID));
+    AliDebug(1, Form("%s: Min TPC PID cut: [%d|%d]\n", GetName(), fMinNClustersTPCPID, nclsTPCPID));
     if(nclsTPCPID >= fMinNClustersTPCPID) SETBIT(survivedCut, kMinNClustersTPCPID);
   }
+  if(TESTBIT(fRequirements, kDriftITS)){
+    //printf("Require drift\n");
+    switch(fDriftITS){
+    case kFirstD:
+      if(TESTBIT(itsPixel, 2)) SETBIT(survivedCut, kDriftITS);
+      break;
+    default: 
+      SETBIT(survivedCut, kDriftITS);
+      break;
+  }
+  }
   if(TESTBIT(fRequirements, kPixelITS)){
     // cut on ITS pixel layers
-    AliDebug(1, "ITS cluster Map: ");
+    AliDebug(1, Form("%s: ITS cluster Map: ", GetName()));
     //PrintBitMap(itsPixel);
     switch(fPixelITS){
       case kFirst:
         AliDebug(2, "First");
+       //printf("First\n");
              if(TESTBIT(itsPixel, 0)) 
                SETBIT(survivedCut, kPixelITS);
         else
@@ -288,6 +369,7 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
                  SETBIT(survivedCut, kPixelITS);
                    break;
       case kSecond: 
+       //printf("Second\n");
         AliDebug(2, "Second");
              if(TESTBIT(itsPixel, 1))
                SETBIT(survivedCut, kPixelITS);
@@ -296,6 +378,7 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
                  SETBIT(survivedCut, kPixelITS);
                    break;
       case kBoth: 
+       //printf("Both\n");
         AliDebug(2, "Both");
              if(TESTBIT(itsPixel,0) && TESTBIT(itsPixel,1))
                      SETBIT(survivedCut, kPixelITS);
@@ -304,6 +387,7 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
                        SETBIT(survivedCut, kPixelITS);
              break;
       case kAny: 
+       //printf("Any\n");
         AliDebug(2, "Any");
              if(TESTBIT(itsPixel, 0) || TESTBIT(itsPixel, 1))
                SETBIT(survivedCut, kPixelITS);
@@ -321,11 +405,17 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
             SETBIT(survivedCut, kPixelITS);
         }
         break;
+      case kNone:
+        // No cut applied, set as survived
+        SETBIT(survivedCut, kPixelITS);
+        break;
       default: 
+        // default, defined as no cut applied 
         AliDebug(2, "None");
+        SETBIT(survivedCut, kPixelITS);
         break;
     }
-    AliDebug(1, Form("Survived Cut: %s\n", TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
+    AliDebug(1, Form("%s: Survived Cut: %s\n", GetName(), TESTBIT(survivedCut, kPixelITS) ? "YES" : "NO"));
   }
 
   if(TESTBIT(fRequirements, kTOFPID)){
@@ -340,24 +430,43 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
 
   if(TESTBIT(fRequirements, kTPCPIDCleanUp)){
       // cut on TPC PID cleanup
-      Int_t fClusterdEdx=GetTPCnclusdEdx(track);
       Bool_t fBitsAboveThreshold=GetTPCCountSharedMapBitsAboveThreshold(track);
-      if((fBitsAboveThreshold==0)&&(fClusterdEdx>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
+      if((fBitsAboveThreshold==0)&&(nclsTPCPID>80)) SETBIT(survivedCut, kTPCPIDCleanUp);
   }
 
   if(TESTBIT(fRequirements, kEMCALmatch)){
     if(track->GetStatus() && AliESDtrack::kEMCALmatch) SETBIT(survivedCut, kEMCALmatch);
   }
 
+  if(TESTBIT(fRequirements, kRejectKinkDaughter)){
+    //printf("test daughter\n");
+    if(!IsKinkDaughter(track)) SETBIT(survivedCut, kRejectKinkDaughter);
+    //else printf("Found kink daughter\n");
+  }
+
+  if(TESTBIT(fRequirements, kRejectKinkMother)){
+    //printf("test mother\n");
+    if(!IsKinkMother(track)) SETBIT(survivedCut, kRejectKinkMother);
+    //else printf("Found kink mother\n");
+  }
+  if(TESTBIT(fRequirements, kTOFsignalDxy)){
+    // cut on TOF matching cluster
+    if((TMath::Abs(tofsignalDx) <= fTOFsignalDx) && (TMath::Abs(tofsignalDz) <= fTOFsignalDz)) SETBIT(survivedCut, kTOFsignalDxy);
+  }
+  if(TESTBIT(fRequirements, kITSpattern)){
+    // cut on ITS pattern (every layer with a working ITS module must have an ITS cluster)
+    if(CheckITSpattern(track)) SETBIT(survivedCut, kITSpattern); 
+  }
+  
   if(fRequirements == survivedCut){
-    //
+    // 
     // Track selected
     //
-    AliDebug(2, "Track Survived cuts\n");
+    AliDebug(2, Form("%s: Track Survived cuts\n", GetName()));
     if(IsQAOn()) FillQAhistosRec(track, kAfterCuts);
     return kTRUE;
   }
-  AliDebug(2, "Track cut");
+  AliDebug(2, Form("%s: Track cut", GetName()));
   if(IsQAOn()) FillCutCorrelation(survivedCut);
   return kFALSE;
 }
@@ -416,7 +525,11 @@ void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
     if(track->GetStatus() && AliESDtrack::kEMCALmatch) hStatusBits->Fill(3);
     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
   }
-  if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(GetTPCnclusdEdx(track));
+  if((htmp = dynamic_cast<TH1F *>(fQAlist->At(7 + when * fgkNQAhistos)))) htmp->Fill(track->GetTPCsignalN());
+
+  if((htmp = dynamic_cast<TH1F *>(fQAlist->At(8 + when * fgkNQAhistos)))) htmp->Fill(GetTRDchi(track));
+
+  if((htmp = dynamic_cast<TH1F *>(fQAlist->At(9 + when * fgkNQAhistos)))) htmp->Fill(GetITSNbOfcls(track));
 }
 
 // //______________________________________________________
@@ -470,7 +583,7 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){
     fQAlist->AddAt(histo1D, 1 + icond * fgkNQAhistos);
     histo1D->GetXaxis()->SetTitle("Impact Parameter");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
-    qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1)), 2 + icond * fgkNQAhistos);
+    qaList->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 100, 0, 1.3)), 2 + icond * fgkNQAhistos);
     fQAlist->AddAt(histo1D, 2 + icond * fgkNQAhistos);
     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
@@ -500,6 +613,15 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){
     fQAlist->AddAt(histo1D, 7 + icond * fgkNQAhistos);
     histo1D->GetXaxis()->SetTitle("Number of TPC clusters for dEdx calculation");
     histo1D->GetYaxis()->SetTitle("counts");
+    qaList->AddAt((histo1D = new TH1F(Form("%s_trdchi2perTracklet%s",GetName(),cutstr[icond].Data()), "chi2 per TRD tracklet", 100, 0, 10)), 8 + icond * fgkNQAhistos);
+    fQAlist->AddAt(histo1D, 8 + icond * fgkNQAhistos);
+    histo1D->GetXaxis()->SetTitle("Chi2 per TRD Tracklet");
+    histo1D->GetYaxis()->SetTitle("Number of Tracks");
+    qaList->AddAt((histo1D = new TH1F(Form("%s_ITScls%s",GetName(),cutstr[icond].Data()), "ITScls", 6, 0, 6)), 9 + icond * fgkNQAhistos);
+    fQAlist->AddAt(histo1D, 9 + icond * fgkNQAhistos);
+    histo1D->GetXaxis()->SetTitle("ITScls");
+    histo1D->GetYaxis()->SetTitle("Number of ITS cls");
+
   }
   // Add cut correlation
   qaList->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2 * fgkNQAhistos);
@@ -535,7 +657,7 @@ Bool_t AliHFEextraCuts::CheckITSstatus(Int_t itsStatus) const {
 }
 
 //______________________________________________________
-Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
+Int_t AliHFEextraCuts::GetITSstatus(const AliVTrack * const track, Int_t layer) const {
        //
        // Check ITS layer status
        //
@@ -543,7 +665,7 @@ Int_t AliHFEextraCuts::GetITSstatus(AliVTrack *track, Int_t layer){
        if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
                Int_t det;
                Float_t xloc, zloc;
-               AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+               const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
                if(esdtrack) esdtrack->GetITSModuleIndexInfo(layer, det, status, xloc, zloc);
        }
        return status;
@@ -556,87 +678,77 @@ Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track)
   // Checks if number of shared bits is above 1
   //
   Int_t nsharebit = 1;
-  TString type = track->IsA()->GetName();
-  if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack){ // coverity
-       TBits shared = esdtrack->GetTPCSharedMap();
-       if(shared.CountBits() >= 2) nsharebit=1;
-        else nsharebit=0;
-    }
-  }
+  const TBits *shared;
+  if(track->IsA() == AliESDtrack::Class())
+    shared = &((static_cast<AliESDtrack *>(track))->GetTPCSharedMap());
+  else if(track->IsA() == AliAODTrack::Class())
+    shared = &((static_cast<AliAODTrack *>(track))->GetTPCSharedMap());
+   else 
+    return 0;
+
+       if(shared->CountBits() >= 2) nsharebit=1;
+  else nsharebit=0;
   return nsharebit;
 }
 
-
-
 //______________________________________________________
 UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
-       //
-       // Get Number of findable clusters in the TPC
-       //
-       Int_t nClusters = 0; // in case no Information available consider all clusters findable
-       TString type = track->IsA()->GetName();
-       if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack){ // coverity
-      if(TESTBIT(fTPCclusterDef, kFoundIter1)){
-                   nClusters = esdtrack->GetTPCNclsIter1();
-        AliDebug(2, ("Using def kFoundIter1"));
-      } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
-        AliDebug(2, ("Using def kCrossedRows"));
-        nClusters = static_cast<UInt_t>(esdtrack->GetTPCClusterInfo(2,1));
-      } else{
-        AliDebug(2, ("Using def kFound"));
-                   nClusters = esdtrack->GetTPCNcls();
-      }
-    }
-  }
-       else if(!type.CompareTo("AliAODTrack")){
-               AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
-    if(aodtrack){
-           const TBits &tpcmap = aodtrack->GetTPCClusterMap();
-           for(UInt_t ibit = 0; ibit < tpcmap.GetNbits(); ibit++)
-                   if(tpcmap.TestBitNumber(ibit)) nClusters++;
-    }
-       }
-       return nClusters;
-}
-
-//______________________________________________________
-UInt_t AliHFEextraCuts::GetTPCnclusdEdx(AliVTrack *track){
   //
-  // Get number of TPC cluster used to calculate dEdx
+  // Get Number of findable clusters in the TPC
   //
-  Int_t nClustersdEdx = 0;
-  TString type = track->IsA()->GetName();
-  if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack){ // coverity
-      nClustersdEdx = esdtrack->GetTPCsignalN();
+  Int_t nClusters = 0; // in case no Information available consider all clusters findable
+  TClass *type = track->IsA();
+  if(TESTBIT(fTPCclusterDef, kFoundIter1)){
+    if(type == AliESDtrack::Class()){
+      AliDebug(2, ("Using def kFoundIter1"));
+      AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+      nClusters = esdtrack->GetTPCNclsIter1();
+    } else {
+      AliDebug(2, ("Number of clusters in iteration 1 not available for AOD tracks"));
     }
+  } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
+    AliDebug(2, ("Using def kCrossedRows"));
+    nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
+  } else if(TESTBIT(fTPCclusterDef, kFound)){
+    AliDebug(2, ("Using def kFound"));
+    nClusters = track->GetTPCNcls();
   }
-  return nClustersdEdx;
+  else {
+    AliDebug(2, ("Using def kFoundAll"));
+    if(type == AliESDtrack::Class()){
+      AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+      const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
+      nClusters = clusterTPC.CountBits();
+    } 
+    else {
+      AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+      const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
+      nClusters = clusterTPC.CountBits();
+    }  
+  }
+  return nClusters;
 }
 
-
-
 //______________________________________________________
 Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
   //
   // Get fraction of shared TPC clusters
   //
   Float_t fracClustersTPCShared = 0.0; 
-  TString type = track->IsA()->GetName();
-  if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack){ // coverity
-      Float_t nClustersTPC = esdtrack->GetTPCclusters(0);
-      Int_t nClustersTPCShared = esdtrack->GetTPCnclsS();
-      if (nClustersTPC!=0) {
-              fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
-      }
-    }
+  Int_t nClustersTPC = track->GetTPCNcls();
+  Int_t nClustersTPCShared = 0;
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    nClustersTPCShared = esdtrack->GetTPCnclsS();
+  } else if(type == AliAODTrack::Class()){
+    AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+    const TBits &shared = aodtrack->GetTPCSharedMap();
+    nClustersTPCShared = shared.CountBits(0) - shared.CountBits(159);
+  }
+  if (nClustersTPC!=0) {
+    fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
   }
   return fracClustersTPCShared;
 }
@@ -648,133 +760,306 @@ Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
   // Only implemented for ESD tracks
   //
   Double_t clusterRatio = 1.; // in case no Information available consider all clusters findable
-  TString type = track->IsA()->GetName();
-  if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack){ // coverity
-      if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
-        AliDebug(2, "Using ratio def kCROverFindable");
-        clusterRatio = esdtrack->GetTPCNclsF() ? esdtrack->GetTPCClusterInfo(2,1)/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // crossed rows/findable
-      } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
-        AliDebug(2, "Using ratio def kFoundOverCR");
-        clusterRatio = esdtrack->GetTPCClusterInfo(2,0); // found/crossed rows
-      } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
-        AliDebug(2, "Using ratio def kFoundOverFindableIter1");
-        clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
-      } else {
-        AliDebug(2, "Using ratio def kFoundOverFindable");
-        clusterRatio = esdtrack->GetTPCNclsF() ? static_cast<Double_t>(esdtrack->GetTPCNcls())/static_cast<Double_t>(esdtrack->GetTPCNclsF()) : 1.; // found/findable
-      }
+  TClass *type = track->IsA();
+  if(TESTBIT(fTPCclusterRatioDef, kCROverFindable)){
+    AliDebug(2, "Using ratio def kCROverFindable");
+    clusterRatio = track->GetTPCNclsF() ? track->GetTPCClusterInfo(2,1)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // crossed rows/findable
+  } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverCR)){
+    AliDebug(2, "Using ratio def kFoundOverCR");
+    clusterRatio = track->GetTPCClusterInfo(2,0); // found/crossed rows
+  } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindableIter1)){
+    if(type == AliESDtrack::Class()){
+      AliDebug(2, "Using ratio def kFoundOverFindableIter1");
+      AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+      clusterRatio = esdtrack->GetTPCNclsFIter1() ? static_cast<Double_t>(esdtrack->GetTPCNclsIter1())/static_cast<Double_t>(esdtrack->GetTPCNclsFIter1()) : 1.; // crossed
+    } else {
+      AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
+      clusterRatio = 1.;
     }
+  } else if(TESTBIT(fTPCclusterRatioDef, kFoundOverFindable)) {
+    AliDebug(2, "Using ratio def kFoundOverFindable");
+    clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // found/findable
+  }
+  else {
+    AliDebug(2, "Using ratio def kFoundAllOverFindable");
+    Int_t allclusters = 0;
+    if(type == AliESDtrack::Class()){
+      AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+      const TBits &clusterTPC = esdtrack->GetTPCClusterMap();
+      allclusters = clusterTPC.CountBits();
+    }
+    else {
+      AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+      const TBits &clusterTPC = aodtrack->GetTPCClusterMap();
+      allclusters = clusterTPC.CountBits();
+    }
+    clusterRatio = track->GetTPCNclsF() ? static_cast<Double_t>(allclusters)/static_cast<Double_t>(track->GetTPCNclsF()) : 1.; // foundall/findable
   }
   return clusterRatio;
+
 }
+//___________________________________________________
+void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
+  //
+  // Get impact parameter
+  //
 
+  const Double_t kBeampiperadius=3.;
+  Double_t dcaD[2]={-999.,-999.},
+           covD[3]={-999.,-999.,-999.};
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    esdtrack->GetImpactParameters(radial, z);
+  }
+  else if(type == AliAODTrack::Class()){
+
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    if(etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+      radial = dcaD[0];
+      z = dcaD[1];
+    }
+    //if(vtxAODSkip) delete vtxAODSkip;
+    if(aodtrack) delete aodtrack;
+  }
+}
 //______________________________________________________
-void AliHFEextraCuts::GetImpactParameters(AliVTrack *track, Float_t &radial, Float_t &z){
-       //
-       // Get impact parameter
-       //
-       TString type = track->IsA()->GetName();
-       if(!type.CompareTo("AliESDtrack")){
-               AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(esdtrack) esdtrack->GetImpactParameters(radial, z);
-       }
-       else if(!type.CompareTo("AliAODTrack")){
-               AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
+  //
+  // Is kink Daughter
+  //
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    if(!esdtrack) return kFALSE;
+    if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
+    else return kFALSE;
+
+  }
+  else if(type == AliAODTrack::Class()){
+    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
     if(aodtrack){
-                 Double_t xyz[3];
-                 aodtrack->XYZAtDCA(xyz);
-                 z = xyz[2];
-                 radial = TMath::Sqrt(xyz[0]*xyz[0] + xyz[1]+xyz[1]);
+      AliAODVertex *aodvertex = aodtrack->GetProdVertex();
+      if(!aodvertex) return kFALSE;
+      if(aodvertex->GetType()==AliAODVertex::kKink) return kTRUE;
+      else return kFALSE;
     }
-       }
+    else return kFALSE;
+  }
+
+  return kFALSE;
+}
+//______________________________________________________
+Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
+  //
+  // Is kink Mother: only for ESD since need to loop over vertices for AOD
+  //
+  //
+
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    if(!esdtrack) return kFALSE;
+    if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
+    else return kFALSE;
+  }
+
+  return kFALSE;
+
+}
+
+//______________________________________________________
+Float_t AliHFEextraCuts::GetTRDchi(AliVTrack *track){
+  //
+  // Get TRDchi2
+  //
+  Int_t ntls(0);
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+      AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+      ntls = esdtrack->GetTRDntracklets();
+      return ntls ? esdtrack->GetTRDchi2()/ntls : -999;
+  }
+  else if(type == AliAODTrack::Class()){
+    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+    if(aodtrack){
+      return  999.;
+    }
+  }
+
+  return 999.;
+
 }
 
+//______________________________________________________
+Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
+  //
+  // Get ITS nb of clusters
+  //
+  TClass *type = track->IsA();
+  if(type == AliESDtrack::Class()){
+    AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
+    return esdtrack->GetITSclusters(0);
+
+  }
+  else if(type == AliAODTrack::Class()){
+    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+    if(aodtrack){
+      return  aodtrack->GetITSNcls();
+    }
+  }
+
+  return 0;
+}
 //______________________________________________________
 void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t &dcaxy, Double_t &dcansigmaxy){
+  //
+  // Get HFE impact parameter (with recalculated primary vertex)
+  //
+  dcaxy=0;
+  dcansigmaxy=0;
+  if(!fEvent){
+    AliDebug(1, "No Input event available\n");
+    return;
+  }
+  TString type = track->IsA()->GetName();
+  const Double_t kBeampiperadius=3.;
+  Double_t dcaD[2]={-999.,-999.},
+           covD[3]={-999.,-999.,-999.};
+
+  if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+    //case of ESD tracks
+    AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+    if(!esdevent) {
+      AliDebug(1, "No esd event available\n");
+      return;
+    }
+
+    //case ESD track: take copy constructor
+    AliESDtrack *esdtrack = NULL;
+    AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
+    if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
+
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+    if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+    }
+
+    if(!vtxESDSkip) return;
+    if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
+      dcaxy = dcaD[0];
+      if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+    }
+    //delete vtxESDSkip;
+    delete esdtrack;
+  }
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+    } 
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    if(etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+      dcaxy = dcaD[0];
+      if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+    }
+    //if(vtxAODSkip) delete vtxAODSkip;
+    if(aodtrack) delete aodtrack;
+  }
+}
+
+//______________________________________________________
+void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
        //
        // Get HFE impact parameter (with recalculated primary vertex)
        //
-       dcaxy=0;
-       dcansigmaxy=0;
   if(!fEvent){
     AliDebug(1, "No Input event available\n");
     return;
   }
   const Double_t kBeampiperadius=3.;
   TString type = track->IsA()->GetName();
-  Double_t dcaD[2]={-999.,-999.};
-  Double_t covD[3]={-999.,-999.,-999.};
-  Float_t dcaF[2]={-999.,-999.};
-  Float_t covF[3]={-999.,-999.,-999.};
   
-  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
-  if(!esdevent) {
-    AliDebug(1, "No esd event available\n");
-    return;
-  }
-  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
-  if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
-   // recalculate primary vertex for peri. and pp
-   AliVertexerTracks vertexer(fEvent->GetMagneticField());
-   vertexer.SetITSMode();
-   vertexer.SetMinClusters(4);
-        Int_t skipped[2];
-   skipped[0] = track->GetID();
-   vertexer.SetSkipTracks(1,skipped);
-   
- //----diamond constraint-----------------------------
-   vertexer.SetConstraintOn();
-   Float_t diamondcovxy[3];
-   esdevent->GetDiamondCovXY(diamondcovxy);
-   Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
-   Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
-   AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
-   vertexer.SetVtxStart(diamond);
-   delete diamond; diamond=NULL;
- //----------------------------------------------------
-    
-   vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
-
-   // Getting the DCA
-   // Propagation always done on a working copy to not disturb the track params of the original track
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-     // Case ESD track: take copy constructor
-     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
-     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack && esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD)){
-    // protection
-    dcaxy = dcaD[0];
-    if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
-    if(!covD[0]) dcansigmaxy = -49.;
-   }
-   delete esdtrack;
-   delete vtxESDSkip;
-  } 
-  else{
-   AliESDtrack *esdtrack = NULL;
-   if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
-    // Case ESD track: take copy constructor
+  if(!TString(track->IsA()->GetName()).CompareTo("AliESDtrack")){
+    //case of ESD tracks
+    AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
+    if(!esdevent) {
+      AliDebug(1, "No esd event available\n");
+      return;
+    }
+
+    //case ESD track: take copy constructor
+    AliESDtrack *esdtrack = NULL;
     AliESDtrack *tmptrack = dynamic_cast<AliESDtrack *>(track);
     if(tmptrack) esdtrack = new AliESDtrack(*tmptrack);
-   } else {
-    // Case AOD track: take different constructor
-    esdtrack = new AliESDtrack(track);
-   }
-   if(esdtrack) esdtrack->GetImpactParameters(dcaF, covF); 
-   dcaxy = dcaF[0];
-   if(covF[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covF[0]);
-   if(!covF[0]) dcansigmaxy = -49.;
-   delete esdtrack;
+
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+    if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+    }
+    if(!vtxESDSkip) return;
+
+    if(esdtrack) esdtrack->PropagateToDCA(vtxESDSkip, fEvent->GetMagneticField(), kBeampiperadius, dcaD, covD);
+
+    //delete vtxESDSkip;
+    delete esdtrack;
   }
-}
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    //Case ESD track: take copy constructor
+    AliAODTrack *aodtrack = NULL;
+    AliAODTrack *tmptrack = dynamic_cast<AliAODTrack *>(track);
+    if(tmptrack) aodtrack = new AliAODTrack(*tmptrack);
 
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+    if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+      vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+    }
+    if(!vtxAODSkip) return;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(aodtrack);
+    etp.PropagateToDCA(vtxAODSkip,aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD);
+
+    //delete vtxAODSkip;
+    delete aodtrack;
+  }
+}
 
 //______________________________________________________
 void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
@@ -808,3 +1093,138 @@ void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maxi
     else maximpactRcut = 9999999999.0;
   }
 }
+//______________________________________________________
+void AliHFEextraCuts::GetTOFsignalDxDz(AliVTrack *track, Double_t &tofsignalDx, Double_t &tofsignalDz){
+  //
+  // TOF matching 
+  //
+  
+  TString type = track->IsA()->GetName();
+  if(!type.CompareTo("AliESDtrack")){
+    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+    if(!esdtrack) return;
+    tofsignalDx = esdtrack->GetTOFsignalDx();
+    tofsignalDz = esdtrack->GetTOFsignalDz();
+  }
+
+}
+
+//______________________________________________________
+Bool_t AliHFEextraCuts::CheckITSpattern(const AliVTrack *const track) const {
+  //
+  // Check if every ITS layer, which has a module which is alive, also
+  // has an ITS cluster
+  //
+  Bool_t patternOK(kTRUE);
+  Int_t status(0);
+  for(Int_t ily = 0; ily < 6; ily++){
+    status = GetITSstatus(track, ily);
+    if(CheckITSstatus(status)){
+      // pixel alive, check whether layer has a cluster
+      if(!TESTBIT(track->GetITSClusterMap(),ily)){
+        // No cluster even though pixel is alive - reject track
+        patternOK = kFALSE;
+        break;
+      }
+    }
+  }
+  return patternOK;
+}
+
+//---------------------------------------------------------------------------
+const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(AliESDEvent *esdevent, AliVTrack *track) {
+  //
+  // This method returns a primary vertex without the daughter tracks of the 
+  // candidate and it recalculates the impact parameters and errors for ESD tracks.
+  // 
+  // The output vertex is created with "new". 
+  //
+
+  const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+
+  AliVertexerTracks vertexer(fEvent->GetMagneticField());
+  vertexer.SetITSMode();
+  vertexer.SetMinClusters(4);
+  Int_t skipped[2];
+  skipped[0] = track->GetID();
+  vertexer.SetSkipTracks(1,skipped);
+
+  //diamond constraint
+  vertexer.SetConstraintOn();
+  Float_t diamondcovxy[3];
+  esdevent->GetDiamondCovXY(diamondcovxy);
+  Double_t pos[3]={esdevent->GetDiamondX(),esdevent->GetDiamondY(),0.};
+  Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+  AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
+  vertexer.SetVtxStart(diamond);
+  delete diamond; diamond=NULL;
+
+  vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
+
+  return vtxESDSkip;
+}
+
+//---------------------------------------------------------------------------
+AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(AliAODEvent *aod, AliVTrack *track) {
+  //
+  // This method returns a primary vertex without the daughter tracks of the 
+  // candidate and it recalculates the impact parameters and errors for AOD tracks.
+  // The output vertex is created with "new". 
+  //
+
+  AliAODVertex *vtxAOD = aod->GetPrimaryVertex();
+  if(!vtxAOD) return 0;
+  TString title=vtxAOD->GetTitle();
+  if(!title.Contains("VertexerTracks")) return 0;
+
+  AliVertexerTracks *vertexer = new AliVertexerTracks(aod->GetMagneticField());
+
+  vertexer->SetITSMode();
+  vertexer->SetMinClusters(3);
+  vertexer->SetConstraintOff();
+
+  if(title.Contains("WithConstraint")) {
+    Float_t diamondcovxy[3];
+    aod->GetDiamondCovXY(diamondcovxy);
+    Double_t pos[3]={aod->GetDiamondX(),aod->GetDiamondY(),0.};
+    Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
+    AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
+    vertexer->SetVtxStart(diamond);
+    delete diamond; diamond=NULL;
+  }
+  Int_t skipped[2]; for(Int_t i=0;i<2;i++) skipped[i]=-1;
+  Int_t id = (Int_t)track->GetID();
+  if(!(id<0)) skipped[0] = id;
+
+  /*// exclude tracks with global constrained parameters
+  Int_t nTracks=aod->GetNumberOfTracks();
+  for(Int_t i=0; i<nTracks; i++){
+    t = aod->GetTrack(i);
+    if(t->TestFilterMask(512)){
+      id = (Int_t)t->GetID();
+      skipped[nTrksToSkip++] = id;
+    }
+  }*/
+
+  vertexer->SetSkipTracks(1,skipped);
+  AliESDVertex *vtxESDNew = vertexer->FindPrimaryVertex(aod);
+
+  delete vertexer; vertexer=NULL;
+
+  if(!vtxESDNew) return 0;
+  if(vtxESDNew->GetNContributors()<=0) {
+    delete vtxESDNew; vtxESDNew=NULL;
+    return 0;
+  }
+
+  // convert to AliAODVertex
+  Double_t pos[3],cov[6],chi2perNDF;
+  vtxESDNew->GetXYZ(pos); // position
+  vtxESDNew->GetCovMatrix(cov); //covariance matrix
+  chi2perNDF = vtxESDNew->GetChi2toNDF();
+  delete vtxESDNew; vtxESDNew=NULL;
+
+  AliAODVertex *vtxAODNew = new AliAODVertex(pos,cov,chi2perNDF);
+
+  return vtxAODNew;
+}