]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGHF/correlationHF/AliDxHFEParticleSelectionEl.cxx
Use configuration script
[u/mrichter/AliRoot.git] / PWGHF / correlationHF / AliDxHFEParticleSelectionEl.cxx
index 55254b14762588920cc0e89a37e417aa975926d7..d93f45eb7a8db3afd119b91144c5fc724000ec7c 100644 (file)
@@ -23,6 +23,9 @@
 /// @brief  D0 selection for D0-HFE correlation
 ///
 #include "AliDxHFEParticleSelectionEl.h"
+#include "AliSelectNonHFE.h"
+#include "AliReducedParticle.h"
+#include "AliESDtrackCuts.h"
 #include "AliVParticle.h"
 #include "AliVEvent.h"
 #include "AliPID.h"
 #include "AliHFEtools.h"
 #include "AliHFEcuts.h"
 #include "AliAODTrack.h"
+#include "AliAODEvent.h"
 #include "AliAnalysisDataSlot.h"
 #include "AliAnalysisDataContainer.h"
 #include "AliAnalysisManager.h"
+#include "AliExternalTrackParam.h"
+#include "AliAODVertex.h"
+#include "AliHFEextraCuts.h"
 #include "AliCFManager.h"
 #include "THnSparse.h"
+#include "AliLog.h"
 #include "TH1F.h"
 #include "TH2F.h"
 #include "TAxis.h"
@@ -58,15 +66,25 @@ AliDxHFEParticleSelectionEl::AliDxHFEParticleSelectionEl(const char* opt)
   : AliDxHFEParticleSelection("electron", opt)
   , fPIDTOFTPC(NULL)
   , fPIDTOF(NULL)
+  , fPIDTPC(NULL)
   , fElectronProperties(NULL)
   , fHistoList(NULL)
   , fCutPidList(NULL)
   , fPIDResponse(NULL)
   , fCuts(NULL)
+  , fSelNHFE(NULL)
+  , fTrackCuts(NULL)
   , fCFM(NULL)
   , fFinalCutStep(kPIDTOFTPC)
   , fInvMassLow(0.1)
-  , fUseInvMassCut(kFALSE)
+  , fUseInvMassCut(kNoInvMass)
+  , fSystem(0)
+  , fTrackNum(0)
+  , fImpactParamCutRadial(3)
+  , fEtaCut(0.8)
+  , fSurvivedCutStep(kNotSelected)
+  , fStoreCutStepInfo(kFALSE)
+
 {
   // constructor
   // 
@@ -98,7 +116,16 @@ AliDxHFEParticleSelectionEl::~AliDxHFEParticleSelectionEl()
   // NOTE: external objects fPID and fCuts are not deleted here
   fPIDTOFTPC=NULL;
   fPIDTOF=NULL;
+  fPIDTPC=NULL;
   fCuts=NULL;
+  if(fSelNHFE){
+    delete fSelNHFE;
+    fSelNHFE=NULL;
+  }
+  if(fTrackCuts){
+    delete fTrackCuts;
+    fTrackCuts=NULL;
+  }
 }
 
 const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
@@ -108,9 +135,10 @@ const char* AliDxHFEParticleSelectionEl::fgkCutBinNames[]={
   "kHFEcutsTOFPID",
   "kHFEcutsTPCPID",
   "kPIDTOF",
+  "kPIDTPC",
   "kPIDTPCTOF",
+  "kINVMASS",
   "Selected e"
-
 };
 
 int AliDxHFEParticleSelectionEl::Init()
@@ -119,6 +147,24 @@ int AliDxHFEParticleSelectionEl::Init()
   // init function
   // 
   int iResult=0;
+  // TODO: think about initializing fSelNHFE and fTrackCuts from the macro
+  //Initialization of invariant mass cut function (AliSelectNonHFE)
+  fSelNHFE= new AliSelectNonHFE("IM","IM");
+  // Cuts for associated track in Inv Mass cut
+  fTrackCuts = new AliESDtrackCuts();
+  fTrackCuts->SetAcceptKinkDaughters(kFALSE);
+  fTrackCuts->SetRequireTPCRefit(kTRUE);
+  fTrackCuts->SetEtaRange(-0.9,0.9);
+  fTrackCuts->SetRequireSigmaToVertex(kTRUE);
+  fTrackCuts->SetMaxChi2PerClusterTPC(4.0);
+  fTrackCuts->SetMinNClustersTPC(80);
+  fTrackCuts->SetPtRange(0.3,1e10);
+
+  fSelNHFE->SetTrackCuts(-3, 3, fTrackCuts);
+  fSelNHFE->SetInvariantMassCut(fInvMassLow);
+  fSelNHFE->SetAlgorithm("KF");
+  fSelNHFE->SetAODanalysis(kTRUE);
+
 
   // Implicit call to InitControlObjects() before setting up CFM and fCuts
   // (if not there)
@@ -148,6 +194,24 @@ int AliDxHFEParticleSelectionEl::Init()
   // TODO: error handling?
   fCuts->Initialize(fCFM);
 
+  //Setting up TPC PID
+  //Add settings for asymmetric cut on nSigma TPC
+  //TODO: have this completely set up from addtask 
+  const int paramSize=4;
+  Double_t params[paramSize];
+  memset(params, 0, sizeof(Double_t)*paramSize);
+  params[0]=-1.;
+  fPIDTPC = new AliHFEpid("hfePidTPC");
+  if(!fPIDTPC->GetNumberOfPIDdetectors()) { 
+    fPIDTPC->AddDetector("TPC",1);
+  }
+  fPIDTPC->ConfigureTPCdefaultCut(NULL, params, 3.);
+  fPIDTPC->InitializePID();
+
+  if(fUseInvMassCut>kNoInvMass)  AliDebug(2,Form("Setting up with invariant mass cut of %f",fInvMassLow));
+  if(fStoreCutStepInfo) AliDebug(2,"will store cut step info");
+  AliDebug(2,Form("Stopping selection after step: %d", fFinalCutStep));
+  
   return 0;
 
 }
@@ -188,6 +252,8 @@ int AliDxHFEParticleSelectionEl::InitControlObjects()
   fHistoList->SetOwner();
 
   // Histogram storing which cuts have been applied to the tracks
+  // TODO: revise the names of the histograms, starting with 'f'
+  // confuses the reader to think of class members
   fHistoList->Add(CreateControlHistogram("fWhichCut","effective cut for a rejected particle", kNCutLabels, fgkCutBinNames));
   fHistoList->Add(CreateControlHistogram("fTPCnClAOD","Nr TPC clusters/track for all tracks", 160, 0., 159.));
   fHistoList->Add(CreateControlHistogram("fTPCnClSingleTrackCuts","Nr TPC clusters/track After SingleTrackCuts", 160, 0., 159.));
@@ -200,12 +266,14 @@ int AliDxHFEParticleSelectionEl::InitControlObjects()
   fHistoList->Add(CreateControl2DHistogram("fdEdx", "dEdx before cuts", dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fdEdxCut", "dEdx after cuts",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fdEdxPidTOF", "dEdx after TOF pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
+  fHistoList->Add(CreateControl2DHistogram("fdEdxPidTPC", "dEdx after TPC pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fdEdxPid", "dEdx after pid",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fdEdxIM", "dEdx after Inv Mass",dEdxBins,"momentum (GeV/c)","dE/dx in TPC (a.u.)"));
 
   // nSigmaTPC vs momentum
   fHistoList->Add(CreateControl2DHistogram("fnSigTPC", "nSigTPC before cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fnSigTPCCut", "nSigmaTPC after cuts",nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
+  fHistoList->Add(CreateControl2DHistogram("fnSigTPCPidTPC", "nSigmaTPC after TPC PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPidTOF", "nSigmaTPC after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fnSigTPCPid", "nSigmaTPC after PID", nSigBins,"momentum (GeV/c)","nSigma in TPC (a.u.)"));
 
@@ -213,6 +281,7 @@ int AliDxHFEParticleSelectionEl::InitControlObjects()
   fHistoList->Add(CreateControl2DHistogram("fnSigTOF", "nSigmaTOF before cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));  
   fHistoList->Add(CreateControl2DHistogram("fnSigTOFCut", "nSigmaTOF after cuts",nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));  
   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPidTOF", "nSigmaTOF after TOF PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
+  fHistoList->Add(CreateControl2DHistogram("fnSigTOFPidTPC", "nSigmaTOF after TPC PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
   fHistoList->Add(CreateControl2DHistogram("fnSigTOFPid", "nSigmaTOF after PID", nSigBins,"momentum (GeV/c)","nSigma in TOF (a.u.)"));
 
   // Invariant mass LS and ULS without cut
@@ -220,14 +289,17 @@ int AliDxHFEParticleSelectionEl::InitControlObjects()
   fHistoList->Add(CreateControlHistogram("fInvMassULS", "Invariant mass ULS", 1000, 0., 0.5));
 
   // Invariant mass LS and ULS without cut, extended x-axis
-  fHistoList->Add(CreateControlHistogram("fInvMassLSExt", "Invariant mass LS", 1000, 0., 5.));
-  fHistoList->Add(CreateControlHistogram("fInvMassULSExt", "Invariant mass ULS", 1000, 0., 5.));
+  fHistoList->Add(CreateControlHistogram("fInvMass2SelLS", "Invariant mass two selected particles LS", 1000, 0., 5.));
+  fHistoList->Add(CreateControlHistogram("fInvMass2SelULS", "Invariant mass two selected particles ULS", 1000, 0., 5.));
 
   // Invariant mass LS and ULS with cut
-  fHistoList->Add(CreateControlHistogram("fInvMassLScut", "Invariant mass LS (cut)", 1000, 0., 0.5));
-  fHistoList->Add(CreateControlHistogram("fInvMassULScut", "Invariant mass ULS (cut)", 1000, 0., 0.5));
+  // TODO: Remove if remove first try at invariant mass
+  fHistoList->Add(CreateControlHistogram("fInvMass2SelLScut", "Invariant mass two selected particles LS (cut)", 1000, 0., 0.5));
+  fHistoList->Add(CreateControlHistogram("fInvMass2SelULScut", "Invariant mass two selected particles ULS (cut)", 1000, 0., 0.5));
 
+  fSelNHFE->SetHistMass((TH1F*)fHistoList->FindObject("fInvMassULS"));
+  fSelNHFE->SetHistMassBack((TH1F*)fHistoList->FindObject("fInvMassLS"));
+  
   AddControlObject(fHistoList);
 
   return AliDxHFEParticleSelection::InitControlObjects();
@@ -274,27 +346,35 @@ TObjArray* AliDxHFEParticleSelectionEl::Select(const AliVEvent* pEvent)
   if (!pEvent) return NULL;
 
   TList* selectedTracks=new TList;
-
+  fSelNHFE->SetPIDresponse(fPIDResponse);
+  fPIDTPC->SetPIDResponse(fPIDResponse);
   TObjArray* finalTracks=new TObjArray;
   if (!finalTracks) return NULL;
   finalTracks->SetOwner(kFALSE); // creating new track objects below
+  TObjArray* finalTracksIM=new TObjArray;
+  if (!finalTracksIM) return NULL;
+  finalTracksIM->SetOwner(kFALSE); // creating new track objects below
 
   int nofTracks=pEvent->GetNumberOfTracks();
+  fTrackNum=0;
   int selectionCode=0;
+  //  int origin
   for (int itrack=0; itrack<nofTracks; itrack++) {
     selectionCode=0;
     AliVParticle* track=pEvent->GetTrack(itrack);
     selectionCode=IsSelected(track,pEvent);
+    fTrackNum++; 
     // If fUseInvMassCut is set to true, only call HistogramParticleProperties
     // (here) for those who didn't pass single track cuts and pid. 
     // if fUseInvMassCut is not set, then call the function for all tracks
     // TODO: Change this when have added invariant mass i thnsparse
-    if(!fUseInvMassCut || selectionCode==0)
-      HistogramParticleProperties(track, selectionCode);
+    //CHANGED FOR OLD INV MASS
     if (selectionCode==0) continue;
+    AliReducedParticle *trackReduced2 =(AliReducedParticle*)CreateParticle(track);
+    finalTracks->Add(trackReduced2);
+    if(fUseInvMassCut!=kInvMassTwoSelected || selectionCode==0)
+      HistogramParticleProperties(trackReduced2, selectionCode);
     selectedTracks->Add(track);
-    //if fUseInvMassCut, fill the array below
-    if(!fUseInvMassCut) finalTracks->Add(CreateParticle(track));
   }
 
   // Calculating invariant mass for electron pairs with same selection criteria
@@ -307,35 +387,43 @@ TObjArray* AliDxHFEParticleSelectionEl::Select(const AliVEvent* pEvent)
       InvMassFilter(selectedTracks, selTrackIndex);
 
       //Only remove electrons if fUseInvMassCut is set
-      if(fUseInvMassCut){
+      if(fUseInvMassCut==kInvMassTwoSelected){
        for(Int_t k=0; k<selectedTracks->GetSize(); k++)
          {
            selectionCode=0; //Reset selectionCode here to be based on invariant mass cut 
+           //On the basis that selectedTracks and finalTracks contain same particles
            AliAODTrack *trackIM=(AliAODTrack*)selectedTracks->At(k);
+           AliReducedParticle *trackReduced=(AliReducedParticle*)finalTracks->At(k);
            if(selTrackIndex[k]==kTRUE)
              {
                ((TH2D*)fHistoList->FindObject("fdEdxIM"))->Fill(trackIM->GetTPCmomentum(), trackIM->GetTPCsignal());
-               finalTracks->Add(CreateParticle(trackIM));
+               finalTracksIM->Add(trackReduced);
                selectionCode=1;
              }
-           HistogramParticleProperties(trackIM, selectionCode);
+           HistogramParticleProperties(trackReduced, selectionCode);
          }
       }
       delete [] selTrackIndex;
       }
       
     }
-  return finalTracks;
+
+  HistogramEventProperties(AliDxHFEParticleSelection::kHistoNrTracksPrEvent,finalTracks->GetEntries());
+  if(fUseInvMassCut!=kInvMassTwoSelected) 
+    return finalTracks; 
+  else
+    return finalTracksIM;
+
 }
 
 int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent* pEvent)
 {
   /// select El candidates
-  // TODO: How to handle MC? would be too much duplicate code if copy entire IsSelected. 
   if(!pEvent){
     AliError("No event information");
     return 0;
-  }
+  }  
+  fSurvivedCutStep=kNotSelected;
 
   AliAODTrack *track=(AliAODTrack*)pEl;
   fCFM->SetRecEventInfo(pEvent);
@@ -345,44 +433,107 @@ int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*
   ((TH2D*)fHistoList->FindObject("fnSigTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
   ((TH1D*)fHistoList->FindObject("fTPCnClAOD"))->Fill(track->GetTPCNcls());
 
+
+  if(fFinalCutStep==kNoCuts){
+    Float_t radial=999;
+    //Float_t z=999;
+    Bool_t acceptTrack=kFALSE;
+
+    //Cut on Eta
+    if(TMath::Abs(track->Eta()) < fEtaCut) {
+      acceptTrack=kTRUE;
+    }
+
+    const Double_t kBeampiperadius=fImpactParamCutRadial;
+    Double_t dcaD[2]={-999.,-999.}, covD[3]={-999.,-999.,-999.};
+
+    AliAODEvent *aodevent = (AliAODEvent*)(pEvent);
+    if(!aodevent) {
+      AliDebug(1, "No aod event available\n");
+      return 0;
+    }
+
+    AliAODVertex *vtxAODSkip  = aodevent->GetPrimaryVertex();
+    if(!vtxAODSkip) return 0;
+    AliExternalTrackParam etp; etp.CopyFromVTrack(track);
+    if(etp.PropagateToDCA(vtxAODSkip, aodevent->GetMagneticField(), kBeampiperadius, dcaD, covD)) {
+      radial = dcaD[0];
+      //z = dcaD[1];
+    }
+
+    // Cut on z  
+    if(acceptTrack && TMath::Abs(radial) < fImpactParamCutRadial){
+      acceptTrack=kTRUE;
+    }
+    
+    if(acceptTrack)
+      fSurvivedCutStep=kNoCuts;
+
+    if(!acceptTrack) return 0;
+
+    // Should only return at this point if you only want to store all tracks
+    AliDebug(2,"Returns after kNoCuts"); 
+    ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); 
+    return 1;
+  }
+
   //--------track cut selection-----------------------
   //Using AliHFECuts:
   // RecKine: ITSTPC cuts  
   if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)){
+    AliDebug(4,"Cut: kStepRecKineITSTPC");
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecKineITSTPC);
     return 0;
   }
-  
+  if(fStoreCutStepInfo) fSurvivedCutStep=kRecKineITSTPC;
+  if(fFinalCutStep==kRecKineITSTPC) {AliDebug(2,"Returns after kStepRecKineITSTPC"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
+
   // RecPrim
   if(!ProcessCutStep(AliHFEcuts::kStepRecPrim, track)) {
+    AliDebug(4,"Cut: kStepRecPrim");
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kRecPrim);
-    return 0;
+    if(!fStoreCutStepInfo) return 0;
+    else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
   }
-  
+  if(fStoreCutStepInfo) fSurvivedCutStep=kRecPrim;
+  if(fFinalCutStep==kRecPrim) {AliDebug(2,"Returns after kRecPrim"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
+
   // HFEcuts: ITS layers cuts
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsITS, track)) {
+    AliDebug(4,"Cut: kStepHFEcutsITS");
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsITS);
-    return 0;
+    if(!fStoreCutStepInfo) return 0;
+    else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
+
   }
+  if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsITS;
+  if(fFinalCutStep==kHFEcutsITS) {AliDebug(2,"Returns after kHFEcutsITS "); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
   
   // HFE cuts: TOF PID and mismatch flag
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTOF, track)) {
+    AliDebug(4,"Cut: kStepHFEcutsTOF");
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTOF);
-    return 0;
+    if(!fStoreCutStepInfo) return 0;
+    else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
   }
+  if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTOF;
+  if(fFinalCutStep==kHFEcutsTOF) {AliDebug(2,"Returns after kHFEcutsTOF"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
   
   // HFE cuts: TPC PID cleanup
   if(!ProcessCutStep(AliHFEcuts::kStepHFEcutsTPC, track)){
+    AliDebug(4,"Cut: kStepHFEcutsTPC");
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kHFEcutsTPC);
-    return 0;
+    if(!fStoreCutStepInfo) return 0;
+    else return 1; //return 1 because it passed cuts above, but not this (no need to go further)
   } 
+  if(fStoreCutStepInfo) fSurvivedCutStep=kHFEcutsTPC;
 
   ((TH2D*)fHistoList->FindObject("fdEdxCut"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
   ((TH2D*)fHistoList->FindObject("fnSigTPCCut"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
   ((TH2D*)fHistoList->FindObject("fnSigTOFCut"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
   ((TH1D*)fHistoList->FindObject("fTPCnClSingleTrackCuts"))->Fill(track->GetTPCNcls());
 
-  if(fFinalCutStep==kHFEcutsTPC) {AliDebug(2,"Returns after track cuts"); return 1;}
+  if(fFinalCutStep==kHFEcutsTPC) {AliDebug(2,"Returns after track cuts"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
 
 
   //--------PID selection-----------------------
@@ -391,8 +542,13 @@ int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*
   hfetrack.SetRecTrack(track);
 
   // TODO: configurable colliding system
-  if(GetSystem()==1) hfetrack.SetPbPb();
+  // TODO: Check problem with PbPb, for now workaround with setting multiplicity to -1
+  if(GetSystem()==1) {
+    hfetrack.SetPbPb();
+    hfetrack.SetCentrality(-1); 
+  }
   else  hfetrack.SetPP();
+  //  hfetrack.SetMulitplicity(ncontribVtx);
 
   // TODO: Put this into a while-loop instead, looping over the number of pid objects in the cut-list?
   // This needs a bit of thinking and finetuning (wrt histogramming)
@@ -400,30 +556,67 @@ int AliDxHFEParticleSelectionEl::IsSelected(AliVParticle* pEl, const AliVEvent*
     ((TH2D*)fHistoList->FindObject("fdEdxPidTOF"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
     ((TH2D*)fHistoList->FindObject("fnSigTPCPidTOF"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
     ((TH2D*)fHistoList->FindObject("fnSigTOFPidTOF"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
-
-    if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF"); return 1;}
+    if(fStoreCutStepInfo){fSurvivedCutStep=kPIDTOF; }
+    if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF");((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
   }
   else{
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOF);
   }
 
+  if(fFinalCutStep==kPIDTOF) {AliDebug(2,"Returns at PIDTOF"); return 0;}
+
+
+  if(fPIDTPC && fPIDTPC->IsSelected(&hfetrack)) {
+    ((TH2D*)fHistoList->FindObject("fdEdxPidTPC"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
+    ((TH2D*)fHistoList->FindObject("fnSigTPCPidTPC"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
+    ((TH2D*)fHistoList->FindObject("fnSigTOFPidTPC"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
+    if(fStoreCutStepInfo){fSurvivedCutStep=kPIDTPC; }
+    if(fFinalCutStep==kPIDTPC) {AliDebug(2,"Returns at PIDTPC"); ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected); return 1;}
+  }
+  else{
+    ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTPC);
+    //return 0; //if rejected by TOF, will also be rejected by TPCTOF
+  }
+  if(fFinalCutStep==kPIDTPC) {AliDebug(2,"Returns at PIDTTPC"); return 0;}
+
   //Combined tof & tpc pid
   if(fPIDTOFTPC && fPIDTOFTPC->IsSelected(&hfetrack)) {
     AliDebug(3,"Inside FilldPhi, electron is selected");
-    ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
     ((TH2D*)fHistoList->FindObject("fdEdxPid"))->Fill(track->GetTPCmomentum(), track->GetTPCsignal());
     ((TH2D*)fHistoList->FindObject("fnSigTPCPid"))->Fill(track->GetTPCmomentum(), fPIDResponse->NumberOfSigmasTPC(track,AliPID::kElectron));
     ((TH2D*)fHistoList->FindObject("fnSigTOFPid"))->Fill(track->P(), fPIDResponse->NumberOfSigmasTOF(track,AliPID::kElectron));
     ((TH1D*)fHistoList->FindObject("fTPCnClTPCTOFPID"))->Fill(track->GetTPCNcls());
+    // Only do this if 
+    if(fStoreCutStepInfo){
+      fSurvivedCutStep=kPIDTOFTPC;
 
-    return 1;
+    }
   }
   else{
     ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kPIDTOFTPC);
-    return 0;
+    AliDebug(4,"Cut: kTPCTOFPID");
+    if(!fStoreCutStepInfo) return 0;
+    else { return 1; }//return 1 because it passed cuts above, but not this (no need to go further)
   }
+  
+  if(fUseInvMassCut==kInvMassSingleSelected)
+    {
+      AliDebug(4,"invmass check");
+      fSelNHFE->FindNonHFE(fTrackNum, track, const_cast<AliVEvent*>(pEvent));
+      if(fSelNHFE->IsLS() || fSelNHFE->IsULS())
+       {
+         //Not selected
+         ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kINVMASS);
+         AliDebug(4,"Cut: Invmass");
+         if(!fStoreCutStepInfo) return 0;
+       }
+      else
+       if(fStoreCutStepInfo) fSurvivedCutStep=kINVMASS;
 
-
+    }
+  ((TH1D*)fHistoList->FindObject("fWhichCut"))->Fill(kSelected);
+  return 1;
+  
 }
 
 void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
@@ -473,13 +666,18 @@ void AliDxHFEParticleSelectionEl::SetCuts(TObject* cuts, int level)
        if(iii==2){
          fPIDTOFTPC=dynamic_cast<AliHFEpid*>(obj);
          if (!fPIDTOFTPC) 
-           AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
+           AliError(Form("(TOFTPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
        }
        if(iii==3){ 
          fPIDTOF=dynamic_cast<AliHFEpid*>(obj);
          if (!fPIDTOF) 
-           AliError(Form("cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
+           AliError(Form("(TOF) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
        }
+       /*if(iii=4){
+         fPIDTPC=dynamic_cast<AliHFEpid*>(obj);
+         if (!fPIDTPC) 
+           AliError(Form("(TPC) cuts object is not of required type AliHFEpid but %s", cuts->ClassName()));
+           }*/
       }
       
     }
@@ -512,24 +710,53 @@ int AliDxHFEParticleSelectionEl::ParseArguments(const char* arguments)
     if(argument.BeginsWith("elreco=")){
       argument.ReplaceAll("elreco=", "");
       int selectionStep=kPIDTOFTPC; //Default
+      if(argument.CompareTo("alltracks")==0) selectionStep=kNoCuts;      
+      if(argument.CompareTo("afterreckineitstpc")==0) selectionStep=kRecKineITSTPC;
+      if(argument.CompareTo("afterrecprim")==0) selectionStep=kRecPrim;
+      if(argument.CompareTo("afterhfeits")==0) selectionStep=kHFEcutsITS;
+      if(argument.CompareTo("afterhfetof")==0) selectionStep=kHFEcutsTOF;
+      if(argument.CompareTo("afterhfetpc")==0) selectionStep=kHFEcutsTPC;
       if(argument.CompareTo("aftertrackcuts")==0) selectionStep=kHFEcutsTPC;
       if(argument.CompareTo("aftertofpid")==0) selectionStep=kPIDTOF;
+      if(argument.CompareTo("aftertpcpid")==0) selectionStep=kPIDTPC;
       if(argument.CompareTo("afterfullpid")==0) selectionStep=kPIDTOFTPC;
       SetFinalCutStep(selectionStep); 
       continue;   
     }
+    if(argument.BeginsWith("twoselectedinvmasscut")){
+      fUseInvMassCut=kInvMassTwoSelected;
+      AliInfo("Using Invariant mass cut for two selected particles");
+      continue;   
+    }
     if(argument.BeginsWith("useinvmasscut")){
-      fUseInvMassCut=kTRUE;
-      AliInfo("Using Invariant mass cut");
+      fUseInvMassCut=kInvMassSingleSelected;
+      AliInfo("Using Invariant mass cut for single selected particle and looser cuts on partner");
       continue;   
     }
     if(argument.BeginsWith("invmasscut=")){
       argument.ReplaceAll("invmasscut=", "");
-      fInvMassLow=argument.Atoi();
-      fUseInvMassCut=kTRUE;
+      fInvMassLow=argument.Atof();
+      AliInfo(Form("Using invariant mass cut: %f",fInvMassLow));
+      //      fUseInvMassCut=kInvMassSingleSelected;
       continue;   
     }
-
+    if(argument.BeginsWith("impactparamcut=")){
+      argument.ReplaceAll("impactparamcut=", "");
+      fImpactParamCutRadial=argument.Atof();
+      AliInfo(Form("Using impact parameter cut: %f",fImpactParamCutRadial));
+      continue;   
+    }
+    if(argument.BeginsWith("etacut=")){
+      argument.ReplaceAll("etacut=", "");
+      fEtaCut=argument.Atof();
+      AliInfo(Form("Using Eta cut: %f",fEtaCut));
+      continue;   
+    }
+    if(argument.BeginsWith("storelastcutstep")){
+      AliInfo("Stores the last cut step");
+      SetStoreLastCutStep(kTRUE);
+      continue;
+    }
     // forwarding of single argument works, unless key-option pairs separated
     // by blanks are introduced
     AliDxHFEParticleSelection::ParseArguments(argument);
@@ -578,23 +805,23 @@ void AliDxHFEParticleSelectionEl::InvMassFilter(TList *elList, Bool_t *selIndx)
              recg.GetMass(mass,width);
 
              if(fFlagLS) {
-               if(mass<0.5)
-                 ((TH1D*)fHistoList->FindObject("fInvMassLS"))->Fill(mass);
-               ((TH1D*)fHistoList->FindObject("fInvMassLSExt"))->Fill(mass);
+               //              if(mass<0.5)
+               // ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
+               ((TH1D*)fHistoList->FindObject("fInvMass2SelLS"))->Fill(mass);
                if(mass>fInvMassLow)
                  {
-                   ((TH1D*)fHistoList->FindObject("fInvMassLScut"))->Fill(mass);
+                   ((TH1D*)fHistoList->FindObject("fInvMass2SelLScut"))->Fill(mass);
                    selIndx[i]=kTRUE;
                    selIndx[j]=kTRUE;
                  }
              }
              if(fFlagULS) {
-               if(mass<0.5)
-                 ((TH1D*)fHistoList->FindObject("fInvMassULS"))->Fill(mass);
-               ((TH1D*)fHistoList->FindObject("fInvMassULSExt"))->Fill(mass);
+               //if(mass<0.5)
+               //  ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);
+               ((TH1D*)fHistoList->FindObject("fInvMass2SelULS"))->Fill(mass);
                if(mass>fInvMassLow)
                  {
-                   ((TH1D*)fHistoList->FindObject("fInvMassULScut"))->Fill(mass);
+                   ((TH1D*)fHistoList->FindObject("fInvMass2SelULScut"))->Fill(mass);
                    selIndx[i]=kTRUE;
                    selIndx[j]=kTRUE;
                  }