]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/hfe/AliHFEextraCuts.cxx
Place the config and root file at the right place
[u/mrichter/AliRoot.git] / PWGHF / hfe / AliHFEextraCuts.cxx
index 81664352529933813348e80c8a14898004f43af4..4d3aac20f77fdfa4773a26cba7571cc713673e9a 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):
@@ -60,13 +63,21 @@ 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),
+  fOppSideIPcut(kFALSE),
+  fTOFsignalDx(1.0),
+  fTOFsignalDz(1.0),
+  fMagField(-10000),
+  fAODFilterBit(0),
+  fListKinkMothers(1000),
+  fNumberKinkMothers(0),
   fCheck(kFALSE),
   fQAlist(0x0) ,
   fDebugLevel(0)
@@ -89,13 +100,21 @@ 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),
+  fOppSideIPcut(c.fOppSideIPcut),
+  fTOFsignalDx(c.fTOFsignalDx),
+  fTOFsignalDz(c.fTOFsignalDz),
+  fMagField(c.fMagField),
+  fAODFilterBit(c.fAODFilterBit),
+  fListKinkMothers(c.fListKinkMothers),
+  fNumberKinkMothers(c.fNumberKinkMothers),
   fCheck(c.fCheck),
   fQAlist(0x0),
   fDebugLevel(0)
@@ -127,13 +146,21 @@ 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;
+    fOppSideIPcut = c.fOppSideIPcut;
+    fTOFsignalDx = c.fTOFsignalDx;
+    fTOFsignalDz = c.fTOFsignalDz;
+    fMagField = c.fMagField;
+    fAODFilterBit = c.fAODFilterBit;
+    fListKinkMothers = c.fListKinkMothers;
+    fNumberKinkMothers = c.fNumberKinkMothers;
     fCheck = c.fCheck;
     fDebugLevel = c.fDebugLevel;
     memcpy(fImpactParamCut, c.fImpactParamCut, sizeof(Float_t) * 4);
@@ -170,7 +197,21 @@ void AliHFEextraCuts::SetRecEventInfo(const TObject *event){
     return ;
   }
   fEvent = (AliVEvent*) event;
-
+  
+  AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+  if(aodevent){
+    // Initialize lookup for kink mother rejection
+    if(aodevent->GetNumberOfVertices() > fListKinkMothers.GetSize()) fListKinkMothers.Set(aodevent->GetNumberOfVertices());
+    fNumberKinkMothers = 0;
+    for(Int_t ivtx = 0; ivtx < aodevent->GetNumberOfVertices(); ivtx++){
+      AliAODVertex *aodvtx = aodevent->GetVertex(ivtx);
+      if(aodvtx->GetType() != AliAODVertex::kKink) continue;
+      AliAODTrack *mother = (AliAODTrack *) aodvtx->GetParent();
+      if(!mother) continue;
+      Int_t idmother = mother->GetID();
+      fListKinkMothers[fNumberKinkMothers++] = idmother;
+    }
+  }
 }
 
 //______________________________________________________
@@ -194,16 +235,17 @@ 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; 
   GetImpactParameters(track, impactR, impactZ);
-  if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
+  if(TESTBIT(fRequirements, kMinHFEImpactParamR) || TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)||TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
     // Protection for PbPb
     GetHFEImpactParameterCuts(track, hfeimpactRcut, hfeimpactnsigmaRcut);
     GetHFEImpactParameters(track, hfeimpactR, hfeimpactnsigmaR);
@@ -216,11 +258,17 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *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);    
@@ -250,29 +298,38 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
     // cut on min. HFE Impact Parameter in Radial direction
     if(TMath::Abs(hfeimpactR) >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamR);
   }
+  if(TESTBIT(fRequirements, kMinHFEImpactParamRcharge)){
+    // cut on min. HFE Impact Parameter in Radial direction, multiplied by particle charge
+    Double_t charge = (Double_t)track->Charge();
+    
+    if(fMagField < 0)charge *= -1.;//the IP distribution side to be chosen depends on magnetic field polarization
+    if(fOppSideIPcut)charge*=-1.;//in case we choose to select electrons from the side of the photon peak
+    Double_t hfeimpactRcharge = hfeimpactR*charge;//switch selected side of the peak
+    if(hfeimpactRcharge >= hfeimpactRcut) SETBIT(survivedCut, kMinHFEImpactParamRcharge);
+  }
   if(TESTBIT(fRequirements, kMinHFEImpactParamNsigmaR)){
     // cut on max. HFE Impact Parameter n sigma in Radial direction
-    if(fAbsHFEImpactParamNsigmaR) {
-//      if((TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) && (TMath::Abs(hfeimpactnsigmaR) < 8)) { //mj debug
-      if(TMath::Abs(hfeimpactnsigmaR) >= hfeimpactnsigmaRcut) {
+    // 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("0: hfeimpactnsigmaR= %lf hfeimpactnsigmaRcut= %lf\n",hfeimpactnsigmaR,hfeimpactnsigmaRcut);
+        //printf("1: 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);
-        }
+    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
@@ -280,41 +337,61 @@ 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("Min ITS clusters: [%d|%d]\n", 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 TPC tracklets
     //printf(Form("Min TPC cut: [%d|%d]\n", fMinNClustersTPC, nclsTPC));
-    AliDebug(1, 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:
@@ -373,7 +450,7 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
         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)){
@@ -406,17 +483,32 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
     //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(TESTBIT(fRequirements, kAODFilterBit)){
+    AliAODTrack *aodtrack = dynamic_cast<AliAODTrack *>(track);
+    if(aodtrack){
+      Int_t aodfilter = 1 << fAODFilterBit;
+      if(aodtrack->TestFilterBit(aodfilter)) SETBIT(survivedCut, kAODFilterBit);
+    }
   }
   
   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;
 }
@@ -476,6 +568,10 @@ void AliHFEextraCuts::FillQAhistosRec(AliVTrack *track, UInt_t when){
     if(GetTPCCountSharedMapBitsAboveThreshold(track)==0) hStatusBits->Fill(4);
   }
   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));
 }
 
 // //______________________________________________________
@@ -529,7 +625,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");
@@ -559,6 +655,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);
@@ -594,7 +699,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
        //
@@ -602,7 +707,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;
@@ -631,13 +736,13 @@ Bool_t AliHFEextraCuts::GetTPCCountSharedMapBitsAboveThreshold(AliVTrack *track)
 
 //______________________________________________________
 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
+  //
+  // Get Number of findable clusters in the TPC
+  //
+  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()){
+    if(type == AliESDtrack::Class()){
       AliDebug(2, ("Using def kFoundIter1"));
       AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
       nClusters = esdtrack->GetTPCNclsIter1();
@@ -647,11 +752,24 @@ UInt_t AliHFEextraCuts::GetTPCncls(AliVTrack *track){
   } else if(TESTBIT(fTPCclusterDef, kCrossedRows)){
     AliDebug(2, ("Using def kCrossedRows"));
     nClusters = static_cast<UInt_t>(track->GetTPCClusterInfo(2,1));
-  } else{
+  } else if(TESTBIT(fTPCclusterDef, kFound)){
     AliDebug(2, ("Using def kFound"));
-         nClusters = track->GetTPCNcls();
+    nClusters = track->GetTPCNcls();
   }
-       return nClusters;
+  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;
 }
 
 //______________________________________________________
@@ -666,13 +784,13 @@ Float_t AliHFEextraCuts::GetTPCsharedClustersRatio(AliVTrack *track){
   if(type == AliESDtrack::Class()){
     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
     nClustersTPCShared = esdtrack->GetTPCnclsS();
-  } else if(type == AliESDtrack::Class()){
+  } 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);
+    fracClustersTPCShared = Float_t(nClustersTPCShared)/Float_t(nClustersTPC);
   }
   return fracClustersTPCShared;
 }
@@ -700,31 +818,65 @@ Double_t AliHFEextraCuts::GetTPCclusterRatio(AliVTrack *track){
       AliDebug(2, "Cluster ratio after iteration 1 not possible for AOD tracks");
       clusterRatio = 1.;
     }
-  } else {
+  } 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()){
-    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]);
+
+    //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;
   }
 }
 //______________________________________________________
@@ -736,8 +888,8 @@ Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
   if(type == AliESDtrack::Class()){
     AliESDtrack *esdtrack = static_cast<AliESDtrack *>(track);
     if(!esdtrack) return kFALSE;
-    if(esdtrack->GetKinkIndex(0)<=0) return kFALSE;
-    else return kTRUE;
+    if(esdtrack->GetKinkIndex(0)>0) return kTRUE;
+    else return kFALSE;
 
   }
   else if(type == AliAODTrack::Class()){
@@ -756,8 +908,7 @@ Bool_t AliHFEextraCuts::IsKinkDaughter(AliVTrack *track){
 //______________________________________________________
 Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
   //
-  // Is kink Mother: only for ESD since need to loop over vertices for AOD
-  //
+  // Is kink Mother 
   //
 
   TClass *type = track->IsA();
@@ -766,11 +917,40 @@ Bool_t AliHFEextraCuts::IsKinkMother(AliVTrack *track){
     if(!esdtrack) return kFALSE;
     if(esdtrack->GetKinkIndex(0)!=0) return kTRUE;
     else return kFALSE;
+  } else if(type == AliAODTrack::Class()){
+    AliAODTrack *aodtrack = static_cast<AliAODTrack *>(track);
+    for(int ikink = 0; ikink < fNumberKinkMothers; ikink++){
+      if(fListKinkMothers[ikink] == aodtrack->GetID()) return kTRUE;
+    }
   }
 
   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){
   //
@@ -792,91 +972,88 @@ Int_t AliHFEextraCuts::GetITSNbOfcls(AliVTrack *track){
   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;
+void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const 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;
   }
-  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 Double_t kBeampiperadius=3.;
+  Double_t dcaD[2]={-999.,-999.},
+           covD[3]={-999.,-999.,-999.};
+  Bool_t isRecalcVertex(kFALSE);
+
+  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;
+    }
+
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+
+    //case ESD track: take copy constructor
+    const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
+    if(tmptrack){
+
+      if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+        vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+        isRecalcVertex = kTRUE;
+      }
+
+      if(vtxESDSkip){
+        AliESDtrack esdtrack(*tmptrack);
+        fMagField = fEvent->GetMagneticField();
+        if(esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD)){
+          dcaxy = dcaD[0];
+          if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+        }
+        if(isRecalcVertex) delete vtxESDSkip;
+      }
+    }
   }
-  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
-    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;
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+
+    //Case ESD track: take copy constructor
+    const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
+    if(tmptrack){ 
+
+      if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+        vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+        isRecalcVertex = kTRUE;
+      } 
+      if(vtxAODSkip){
+        AliAODTrack aodtrack(*tmptrack);
+        AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
+        fMagField = aodevent->GetMagneticField();
+        if(etp.PropagateToDCA(vtxAODSkip,fMagField, kBeampiperadius, dcaD, covD)) {
+          dcaxy = dcaD[0];
+          if(covD[0]) dcansigmaxy = dcaxy/TMath::Sqrt(covD[0]);
+        }
+        if(isRecalcVertex) delete vtxAODSkip;
+      }
+    }
   }
 }
 
 //______________________________________________________
-void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2], Double_t covD[3]){
+void AliHFEextraCuts::GetHFEImpactParameters(const AliVTrack * const track, Double_t dcaD[2], Double_t covD[3]){
        //
        // Get HFE impact parameter (with recalculated primary vertex)
        //
@@ -886,101 +1063,219 @@ void AliHFEextraCuts::GetHFEImpactParameters(AliVTrack *track, Double_t dcaD[2],
   }
   const Double_t kBeampiperadius=3.;
   TString type = track->IsA()->GetName();
-  Float_t dcaF[2]={-999.,-999.};
-  Float_t covF[3]={-999.,-999.,-999.};
+  Bool_t isRecalcVertex(kFALSE);
   
-  AliESDEvent *esdevent = dynamic_cast<AliESDEvent *>(fEvent);
-  if(!esdevent) {
-    AliDebug(1, "No esd event available\n");
-    return;
+  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;
+    }
+
+    // Check whether primary vertex is available
+    const AliVVertex *vtxESDSkip = esdevent->GetPrimaryVertex();
+    if(!vtxESDSkip) return;
+
+    const AliESDtrack *tmptrack = dynamic_cast<const AliESDtrack *>(track);
+    if(tmptrack){
+      //case ESD track: take copy constructor
+
+      if( vtxESDSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+        vtxESDSkip = RemoveDaughtersFromPrimaryVtx(esdevent, track);
+        isRecalcVertex = kTRUE;
+      }
+      if(vtxESDSkip){
+        AliESDtrack esdtrack(*tmptrack);
+        fMagField = fEvent->GetMagneticField();
+        esdtrack.PropagateToDCA(vtxESDSkip, fMagField, kBeampiperadius, dcaD, covD);
+
+        if(isRecalcVertex) delete vtxESDSkip;
+      }
+    }
   }
-  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);
-   delete esdtrack;
-   delete vtxESDSkip;
-  } 
-  else{
-   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->GetImpactParameters(dcaF, covF); 
-   dcaD[0] = dcaF[0];
-   dcaD[1] = dcaF[1];
-   covD[0] = covF[0];
-   covD[1] = covF[1];
-   covD[2] = covF[2];
-   delete esdtrack;
+  else {
+    //case of AOD tracks
+    AliAODEvent *aodevent = dynamic_cast<AliAODEvent *>(fEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return;
+    }
+
+    // Check whether primary vertex is available
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return;
+
+    //Case ESD track: take copy constructor
+    const AliAODTrack *tmptrack = dynamic_cast<const AliAODTrack *>(track);
+    if(tmptrack){
+
+      if(vtxAODSkip->GetNContributors() < 30){ // if vertex contributor is smaller than 30, recalculate the primary vertex
+        vtxAODSkip = RemoveDaughtersFromPrimaryVtx(aodevent, track);
+        isRecalcVertex = kTRUE;
+      }
+      if(vtxAODSkip){
+        AliAODTrack aodtrack(*tmptrack);
+        AliExternalTrackParam etp; etp.CopyFromVTrack(&aodtrack);
+        fMagField = aodevent->GetMagneticField();
+        etp.PropagateToDCA(vtxAODSkip, fMagField, kBeampiperadius, dcaD, covD);
+
+        if(isRecalcVertex) delete vtxAODSkip;
+      }
+    }
   }
 }
 
 //______________________________________________________
-void AliHFEextraCuts::GetHFEImpactParameterCuts(AliVTrack *track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
+void AliHFEextraCuts::GetHFEImpactParameterCuts(const AliVTrack * const track, Double_t &hfeimpactRcut, Double_t &hfeimpactnsigmaRcut){
        //
        // Get HFE impact parameter cut(pt dependent)
        //
   
-  TString type = track->IsA()->GetName();
-  if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
-    if(!esdtrack) return;
-    Double_t pt = esdtrack->Pt();      
-    hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
-    hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
-  }
+  Double_t pt = track->Pt();   
+  hfeimpactRcut = fIPcutParam[0]+fIPcutParam[1]*exp(fIPcutParam[2]*pt);  // abs R cut
+  hfeimpactnsigmaRcut = fIPcutParam[3];                                  // sigma cut
 }
 //______________________________________________________
-void AliHFEextraCuts::GetMaxImpactParameterCutR(AliVTrack *track, Double_t &maximpactRcut){
+void AliHFEextraCuts::GetMaxImpactParameterCutR(const AliVTrack * const track, Double_t &maximpactRcut){
        //
        // Get max impact parameter cut r (pt dependent)
        //
   
+  Double_t pt = track->Pt();   
+  if(pt > 0.15) {
+    maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
+  }
+  else maximpactRcut = 9999999999.0;
+}
+//______________________________________________________
+void AliHFEextraCuts::GetTOFsignalDxDz(const AliVTrack * const track, Double_t &tofsignalDx, Double_t &tofsignalDz){
+  //
+  // TOF matching 
+  //
+  
   TString type = track->IsA()->GetName();
   if(!type.CompareTo("AliESDtrack")){
-    AliESDtrack *esdtrack = dynamic_cast<AliESDtrack *>(track);
+    const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
     if(!esdtrack) return;
-    Double_t pt = esdtrack->Pt();      
-    if(pt > 0.15) {
-      maximpactRcut = 0.0182 + 0.035/TMath::Power(pt,1.01);  // abs R cut
+    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;
+      }
     }
-    else maximpactRcut = 9999999999.0;
   }
+  return patternOK;
+}
+
+//---------------------------------------------------------------------------
+const AliVVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliESDEvent * const esdevent, const AliVTrack * const 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(pos,cov,1.,1);
+  vertexer.SetVtxStart(&diamond);
+
+  const AliVVertex *vtxESDSkip = vertexer.FindPrimaryVertex(fEvent);
+
+  return vtxESDSkip;
+}
+
+//---------------------------------------------------------------------------
+AliAODVertex* AliHFEextraCuts::RemoveDaughtersFromPrimaryVtx(const AliAODEvent * const aod, const AliVTrack * const 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(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(pos,cov,1.,1);
+    vertexer.SetVtxStart(&diamond);
+  }
+  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);
+
+  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;
 }