]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG/FLOW/Tasks/AliAnalysisTaskPhiFlow.cxx
Merge branch 'master' into dev
[u/mrichter/AliRoot.git] / PWG / FLOW / Tasks / AliAnalysisTaskPhiFlow.cxx
index e9ade369986f31a5aa91a6d1587d35175b5e9f25..4d172696ff0c448425dfc0708675c93a5b65edd7 100644 (file)
 // new in this version (use with caution): vzero event plane, event mixing
 
 #include "TChain.h"
-#include "TTree.h"
 #include "TH1.h"
 #include "TH1F.h"
 #include "TH2F.h"
-#include "TCanvas.h"
 #include "TMath.h"
 #include "TObjArray.h"
-#include "AliAnalysisTask.h"
 #include "AliAnalysisTaskSE.h"
 #include "AliAnalysisManager.h"
-#include "AliESDEvent.h"
 #include "AliAODEvent.h"
-#include "AliESDInputHandler.h"
 #include "AliAODInputHandler.h"
 #include "AliCentrality.h"
-#include "AliVEvent.h"
 #include "AliAnalysisTaskPhiFlow.h"
 #include "AliFlowBayesianPID.h"
-#include "AliStack.h"
-#include "AliESDtrackCuts.h"
+#include "AliPIDCombined.h"
 #include "AliMCEvent.h"
 #include "TProfile.h"
 #include "AliFlowCandidateTrack.h"
@@ -49,9 +42,7 @@
 #include "AliFlowCommonConstants.h"
 #include "AliFlowEvent.h"
 #include "TVector3.h"
-#include "AliESDVZERO.h"
 #include "AliAODVZERO.h"
-#include "AliPID.h"
 #include "AliPIDResponse.h"
 #include "AliAODMCParticle.h"
 #include "AliAnalysisTaskVnV0.h"
@@ -66,7 +57,7 @@ ClassImp(AliAnalysisTaskPhiFlow)
 ClassImp(AliPhiMesonHelperTrack)
 
 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
-   fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fAODAnalysis(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fCentrality(999), fVertex(999), fESD(0), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0),fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0)
+   fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0),fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
 {
    // Default constructor
    for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
@@ -83,7 +74,7 @@ AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow() : AliAnalysisTaskSE(),
 }
 //_____________________________________________________________________________
 AliAnalysisTaskPhiFlow::AliAnalysisTaskPhiFlow(const char *name) : AliAnalysisTaskSE(name),
-   fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fAODAnalysis(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fOldTrackParam(0), fRequireTPCStandAlone(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fCentrality(999), fVertex(999), fESD(0), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethod(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0), fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0), fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0)
+   fDebug(0), fIsMC(0), fEventMixing(0), fTypeMixing(0), fQA(0), fV0(0), fMassBins(1), fMinMass(-1.), fMaxMass(0.), fCutsRP(NULL), fNullCuts(0), fPIDResponse(0), fFlowEvent(0), fBayesianResponse(0), fCandidates(0), fCandidateEtaPtCut(0), fCandidateMinEta(0), fCandidateMaxEta(0), fCandidateMinPt(0), fCandidateMaxPt(0), fNPtBins(18), fCentrality(999), fVertex(999), fAOD(0), fPoolManager(0), fOutputList(0), fEventStats(0), fCentralityPass(0), fCentralityNoPass(0), fNOPID(0), fPIDk(0), fNOPIDTOF(0), fPIDTOF(0), fPtP(0), fPtN(0), fPtKP(0), fPtKN(0), fMultCorAfterCuts(0), fMultvsCentr(0), fCentralityMin(0), fCentralityMax(100), fkCentralityMethodA(0), fkCentralityMethodB(0), fCentralityCut2010(0), fCentralityCut2011(0), fPOICuts(0), fVertexRange(0), fPhi(0), fPt(0), fEta(0), fVZEROA(0), fVZEROC(0), fTPCM(0)/*, fDeltaDipAngle(0), fDeltaDipPt(0), fApplyDeltaDipCut(0)*/, fDCAAll(0), fDCAXYQA(0), fDCAZQA(0), fDCAPrim(0), fDCASecondaryWeak(0), fDCAMaterial(0), fSubEventDPhiv2(0), fSkipEventSelection(0), fUsePidResponse(0), fPIDCombined(0)
 {
    // Constructor
    for(Int_t i(0); i < 7; i++) fPIDConfig[i] = 0.;
@@ -112,6 +103,7 @@ AliAnalysisTaskPhiFlow::~AliAnalysisTaskPhiFlow()
    if (fFlowEvent) delete fFlowEvent;
    if (fBayesianResponse) delete fBayesianResponse;
    if (fEventMixing) delete fPoolManager;
+   if (fPIDCombined) delete fPIDCombined;
    if (fDebug > 0) cout << " === Deleted instance of AliAnalysisTaskPhiFlow === " << endl;
 }
 //_____________________________________________________________________________
@@ -215,6 +207,12 @@ void AliAnalysisTaskPhiFlow::AddPhiIdentificationOutputObjects()
        fOutputList->Add(fDCAXYQA);
        fDCAZQA = new TH1F("fDCAZQA", "fDCAZQA", 1000, -5, 5);
        fOutputList->Add(fDCAZQA);
+       if(fCentralityCut2010 || fCentralityCut2011) {
+           fMultCorAfterCuts = new TH2F("fMultCorAfterCuts", "TPC vs Global multiplicity (After cuts); Global multiplicity; TPC multiplicity", 100, 0, 3000, 100, 0, 3000);
+           fOutputList->Add(fMultCorAfterCuts);
+           fMultvsCentr = new TH2F("fMultvsCentr", "Multiplicity vs centrality; centrality; Multiplicity", 9, -0.5, 100.5, 101, 0, 3000);
+           fOutputList->Add(fMultvsCentr);
+       }
    }
    if(fIsMC || fQA) {
        fDCAAll = new TH2F("fDCAAll", "fDCAAll", 1000, 0, 10, 1000, -5, 5);
@@ -249,6 +247,13 @@ void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
    if(fDebug > 0) cout << " *** UserCreateOutputObjects() *** " << endl;
    fNullCuts = new AliFlowTrackCuts("null_cuts");
    fBayesianResponse = new AliFlowBayesianPID();
+      // combined pid
+   fPIDCombined = new AliPIDCombined;
+   fPIDCombined->SetDefaultTPCPriors();
+   fPIDCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+
+   // flag to mc
+   if(fSkipEventSelection || fIsMC) fBayesianResponse->SetMC(kTRUE);
    Double_t t(0);
    for(Int_t i = 0; i < 7; i++) t+=TMath::Abs(fPIDConfig[i]);
    if(t < 0.1) AliFatal("No valid PID procedure recognized -- terminating analysis !!!");
@@ -260,7 +265,7 @@ void AliAnalysisTaskPhiFlow::UserCreateOutputObjects()
    cc->SetNbinsMass(fMassBins);  cc->SetNbinsEta(200);           (fMassBins == 1) ? cc->SetNbinsPt(15) : cc->SetNbinsPt(100); // high pt
    cc->SetMassMin(fMinMass);     cc->SetEtaMin(-5.0);            cc->SetPtMin(0);
    cc->SetMassMax(fMaxMass);     cc->SetEtaMax(+5.0);            (fMassBins == 1) ? cc->SetPtMax(15) : cc->SetPtMax(10); // high pt
-   if (!fOldTrackParam) fBayesianResponse->SetNewTrackParam();
+   fBayesianResponse->SetNewTrackParam();
    AliAnalysisManager *man = AliAnalysisManager::GetAnalysisManager();
    if (man) {
       AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
@@ -318,6 +323,7 @@ template <typename T> Double_t AliAnalysisTaskPhiFlow::InvariantMass(const T* tr
    return TMath::Sqrt((es - (pxs + pys + pzs)));
 }
 //_____________________________________________________________________________
+/*
 template <typename T> Double_t AliAnalysisTaskPhiFlow::DeltaDipAngle(const T* track1, const T* track2) const
 {
    // Calculate the delta dip angle between two particles
@@ -333,6 +339,7 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckDeltaDipAngle(const T*
    if ((TMath::Abs(DeltaDipAngle(track1, track2)) < fDeltaDipAngle) && (PhiPt(track1, track2) < fDeltaDipPt)) return kFALSE;
    return kTRUE;
 }
+*/
 //_____________________________________________________________________________
 template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCandidateEtaPtCut(const T* track1, const T* track2) const
 {
@@ -351,6 +358,7 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::EventCut(T* event)
    // Impose event cuts
    if(fDebug > 0) cout << " *** EventCut() *** " << endl;
    if (!event) return kFALSE;
+   if (fSkipEventSelection) return kTRUE;
    if (!CheckVertex(event)) return kFALSE;
    if (!CheckCentrality(event)) return kFALSE;
    if(fQA) PlotMultiplcities(event);
@@ -380,23 +388,80 @@ template <typename T> Bool_t AliAnalysisTaskPhiFlow::CheckCentrality(T* event)
 {
    // Check if event is within the set centrality range. Falls back to V0 centrality determination if no method is set
    if(fDebug > 0) cout << " *** CheckCentrality() *** " << endl;
-   if (!fkCentralityMethod) AliFatal("No centrality method set! FATAL ERROR!");
-   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethod);
-   if ((fCentrality <= fCentralityMin) || (fCentrality > fCentralityMax)) {
+   if (!fkCentralityMethodA) AliFatal("No centrality method set! FATAL ERROR!");
+   fCentrality = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodA);
+   Double_t cenB(-999);
+   // if a second centrality estimator is requited, set it
+   (fkCentralityMethodB) ? cenB = event->GetCentrality()->GetCentralityPercentile(fkCentralityMethodB) : cenB = fCentrality;
+   if (TMath::Abs(fCentrality-cenB) > 5 || cenB >= 80 || cenB < 0 || fCentrality <= fCentralityMin || fCentrality > fCentralityMax) {
       if(fQA) fCentralityNoPass->Fill(fCentrality) ;
       return kFALSE;
    }
+   const Int_t nGoodTracks = event->GetNumberOfTracks();
+   if(fCentralityCut2010) { // cut on outliers
+      Float_t multTPC(0.); // tpc mult estimate
+      Float_t multGlob(0.); // global multiplicity
+      for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
+          AliAODTrack* trackAOD = event->GetTrack(iTracks);
+          if (!trackAOD) continue;
+          if (!(trackAOD->TestFilterBit(1))) continue;
+          if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
+          multTPC++;
+      }
+      for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
+          AliAODTrack* trackAOD = event->GetTrack(iTracks);
+          if (!trackAOD) continue;
+          if (!(trackAOD->TestFilterBit(16))) continue;
+          if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
+          Double_t b[2] = {-99., -99.};
+          Double_t bCov[3] = {-99., -99., -99.};
+          AliAODTrack copy(*trackAOD);
+          if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
+          if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
+          multGlob++;
+      } //track loop
+ //     printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
+      if(! (multTPC > (-40.3+1.22*multGlob) && multTPC < (32.1+1.59*multGlob))) return kFALSE;
+      if(fQA) {  
+          fMultCorAfterCuts->Fill(multGlob, multTPC);  
+          fMultvsCentr->Fill(fCentrality, multTPC);
+      }
+   }
+
+   if(fCentralityCut2011) { // cut on outliers
+      Float_t multTPC(0.); // tpc mult estimate
+      Float_t multGlob(0.); // global multiplicity
+      for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill tpc mult
+          AliAODTrack* trackAOD = event->GetTrack(iTracks);
+          if (!trackAOD) continue;
+          if (!(trackAOD->TestFilterBit(1))) continue;
+          if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70)  || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.2)) continue;
+          multTPC++;
+      }
+      for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { // fill global mult
+          AliAODTrack* trackAOD = event->GetTrack(iTracks);
+          if (!trackAOD) continue;
+          if (!(trackAOD->TestFilterBit(16))) continue;
+          if ((trackAOD->Pt() < .2) || (trackAOD->Pt() > 5.0) || (TMath::Abs(trackAOD->Eta()) > .8) || (trackAOD->GetTPCNcls() < 70) || (trackAOD->GetDetPid()->GetTPCsignal() < 10.0) || (trackAOD->Chi2perNDF() < 0.1)) continue;
+          Double_t b[2] = {-99., -99.};
+          Double_t bCov[3] = {-99., -99., -99.};
+          AliAODTrack copy(*trackAOD);
+          if (!(copy.PropagateToDCA(event->GetPrimaryVertex(), event->GetMagneticField(), 100., b, bCov))) continue;
+          if ((TMath::Abs(b[0]) > 0.3) || (TMath::Abs(b[1]) > 0.3)) continue;
+          multGlob++;
+      } //track loop
+      //printf(" mult TPC %.2f, mult Glob %.2f \n", multTPC, multGlob);
+      if(! (multTPC > (-36.73 + 1.48*multGlob) && multTPC < (62.87 + 1.78*multGlob))) return kFALSE;
+      if(fQA) {  
+          fMultCorAfterCuts->Fill(multGlob, multTPC);  
+          fMultvsCentr->Fill(fCentrality, multTPC);
+      }
+   }
+
    fCentralityPass->Fill(fCentrality);
    return kTRUE;
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliESDEvent* event)
-{
-   // Initialize the Bayesian PID object for ESD analysis
-   if(fDebug > 0) cout << " *** InitializeBayesianPID() *** " << endl;
-   fBayesianResponse->SetDetResponse(event, fCentrality, AliESDpid::kTOF_T0, kTRUE);
-}
-//_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::InitializeBayesianPID(AliAODEvent* event)
 {
    // Initialize the Bayesian PID object for AOD
@@ -427,20 +492,22 @@ Bool_t AliAnalysisTaskPhiFlow::PassesDCACut(AliAODTrack* track) const
     // fDCAConfig[0] >  1 pt dependent dca cut
     if(fDebug > 1) cout << " *** PassesDCACut() *** " << endl;
     if(fIsMC) return kTRUE;
+    AliAODTrack copy(*track);
     if( (fDCAConfig[0] < 0.1) && (fDCAConfig[0] > -0.1) ) {
         if(fQA) {
             Double_t b[2] = { -99., -99.};
             Double_t bCov[3] = { -99., -99., -99.};
-            track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
-            fDCAXYQA->Fill(b[0]);
-            fDCAZQA->Fill(b[1]);
-            fDCAPrim->Fill(track->Pt(), b[0]);
+            if(copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) {
+                fDCAXYQA->Fill(b[0]);
+                fDCAZQA->Fill(b[1]);
+                fDCAPrim->Fill(track->Pt(), b[0]);
+            }
         }
         return kTRUE;
     }
     Double_t b[2] = { -99., -99.};
     Double_t bCov[3] = { -99., -99., -99.};
-    track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+    if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return kFALSE;
     if((!fIsMC)&&fQA) fDCAMaterial->Fill(track->Pt(), b[0]);
     if( (fDCAConfig[0] < -.9) && ( (TMath::Abs(b[0]) > fDCAConfig[1]) || (TMath::Abs(b[1]) > fDCAConfig[2])) ) return kFALSE;
     if(fDCAConfig[0] > .9) {
@@ -457,25 +524,15 @@ Bool_t AliAnalysisTaskPhiFlow::PassesDCACut(AliAODTrack* track) const
     return kTRUE;
 }
 //_____________________________________________________________________________
-Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliESDtrack* track) const
-{
-   // Check if particle is a kaon according to method set in steering macro
-   if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
-   if (fRequireTPCStandAlone && (track->GetStatus()&AliESDtrack::kTPCin) == 0) return kFALSE;
-   if (fPIDConfig[1]!=0 || fPIDConfig[4]!=0) AliFatal("TPC || ITS PID not available in ESD anlaysis -- terminating analysis !!!");
-   if (PassesTPCbayesianCut(track))
-   {
-     if(fQA) fPIDk->Fill(track->P(), track->GetTPCsignal());
-     return kTRUE;
-   }
-   return kFALSE;
-}
-//_____________________________________________________________________________
 Bool_t AliAnalysisTaskPhiFlow::IsKaon(AliAODTrack* track) const
 {
    // Kaon identification routine, based on multiple detectors and approaches
    if(fDebug > 1) cout << " *** IsKaon() *** " << endl;
-   if(fRequireTPCStandAlone && (!track->TestFilterBit(1))) return kFALSE;
+   if(fUsePidResponse) {
+       Double_t prob[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
+       fPIDCombined->ComputeProbabilities(track, fPIDResponse, prob);
+       if(prob[3] > fPIDConfig[6])  return kTRUE;
+   }
    if(!PassesDCACut(track)) return kFALSE;
    if(fQA) {fNOPID->Fill(track->P(), track->GetTPCsignal());fNOPIDTOF->Fill(track->P(), track->GetTOFsignal());}
    if(track->Pt() < fPIDConfig[1]) {
@@ -652,7 +709,6 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
    }
    fAOD = dynamic_cast<AliAODEvent*>(InputEvent()); // check for aod data type
    if (fAOD) {
-      fAODAnalysis = kTRUE;
       if (!EventCut(fAOD)) return; // check for event cuts
       InitializeBayesianPID(fAOD); // init event objects
       SetNullCuts(fAOD);
@@ -690,7 +746,7 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
       }
       for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
          for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
-            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
+//            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
             if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
             PtSelector(0, up[pTracks], un[nTracks]);
             Double_t pt = PhiPt(up[pTracks], un[nTracks]);
@@ -733,14 +789,14 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
       if(!fEventMixing) { // combinatorial background
           for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
              for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
-                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
+//                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
                 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
                 PtSelector(1, up[pTracks], up[nTracks]);
              }
           }
           for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
              for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
-                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
+//                if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
                 if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
                 PtSelector(2, un[nTracks], un[pTracks]);
              }
@@ -750,107 +806,14 @@ void AliAnalysisTaskPhiFlow::UserExec(Option_t *)
       PostData(1, fOutputList);
       PostData(2, fFlowEvent);
    }
-   fESD = dynamic_cast<AliESDEvent*>(InputEvent()); // check for esd data type (obsolete)
-   if (fESD) {
-      if(fV0 || fEventMixing) {
-          if (fDebug > 0 ) cout << " FATAL ERROR: AOD only analysis called on ESD data, skipping input !!! " << endl;
-          return;
-      }
-      fAODAnalysis = kFALSE;
-      if (!EventCut(fESD)) return; // check event cuts
-      InitializeBayesianPID(fESD);
-      SetNullCuts(fESD);
-      PrepareFlowEvent(fESD->GetNumberOfTracks());
-      if(fQA) fEventStats->Fill(0);
-      Int_t unTracks = fESD->GetNumberOfTracks();
-      AliESDtrack* un[unTracks];
-      AliESDtrack* up[unTracks];
-      Int_t unp(0);
-      Int_t unn(0);
-      for (Int_t iTracks = 0; iTracks < unTracks; iTracks++) { // select analysis candidates
-         AliESDtrack* track = fESD->GetTrack(iTracks);
-         if (!PhiTrack(track)) continue;
-         Bool_t charge = kFALSE;
-         if (track->Charge() > 0) {
-            charge = kTRUE;
-            if(fQA) {fEventStats->Fill(1); fPtP->Fill(track->Pt());}
-         }
-         if (track->Charge() < 0) {
-            if(fQA) fEventStats->Fill(2);
-            fPtN->Fill(track->Pt());
-         }
-         if (IsKaon(track)) {
-            if (charge) {
-               up[unp] = track;
-               unp++;
-               if(fQA) {fEventStats->Fill(3);fPtKP->Fill(track->Pt());}
-            }
-            if (!charge) {
-               un[unn] = track;
-               unn++;
-               if(fQA) {fEventStats->Fill(4); fPtKN->Fill(track->Pt());}
-            }
-         }
-      }
-        for (Int_t pTracks = 0; pTracks < unp ; pTracks++) { // perform analysis
-         for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
-            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], un[nTracks]))) continue;
-            if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(up[pTracks], un[nTracks]))) continue;
-            PtSelector(0, up[pTracks], un[nTracks]);
-            Double_t pt = PhiPt(up[pTracks], un[nTracks]);
-            Double_t mass = InvariantMass(up[pTracks], un[nTracks]);
-            TVector3 a(up[pTracks]->Px(), up[pTracks]->Py(), up[pTracks]->Pz());
-            TVector3 b(un[nTracks]->Px(), un[nTracks]->Py(), up[pTracks]->Pz());
-            TVector3 c = a + b;
-            Double_t phi = c.Phi();
-            Double_t eta = c.Eta();
-            Int_t nIDs[2];
-            nIDs[0] = up[pTracks]->GetID();
-            nIDs[1] = un[nTracks]->GetID();
-            MakeTrack(mass, pt, phi, eta, 2, nIDs);
-         }
-      }
-      if (fDebug > 0)  printf("I received %d candidates\n", fCandidates->GetEntriesFast()); // insert candidates in flow event
-      for (int iCand = 0; iCand != fCandidates->GetEntriesFast(); ++iCand) {
-         AliFlowCandidateTrack *cand = dynamic_cast<AliFlowCandidateTrack*>(fCandidates->At(iCand));
-         if (!cand) continue;
-         if (fDebug > 1) printf(" >Checking at candidate %d with %d daughters: mass %f\n", iCand, cand->GetNDaughters(), cand->Mass());
-         for (int iDau = 0; iDau != cand->GetNDaughters(); ++iDau) {
-            if (fDebug > 1) printf("  >Daughter %d with fID %d", iDau, cand->GetIDDaughter(iDau));
-            for (int iRPs = 0; iRPs != fFlowEvent->NumberOfTracks(); ++iRPs) {
-               AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack(iRPs));
-               if (!iRP) continue;
-               if (!iRP->InRPSelection()) continue;
-               if (cand->GetIDDaughter(iDau) == iRP->GetID()) {
-                  if (fDebug>1) printf(" was in RP set");
-                  iRP->SetForRPSelection(kFALSE);
-                  fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
-               }
-            }
-            if (fDebug>1) printf("\n");
-         }
-         cand->SetForPOISelection(kTRUE);
-         fFlowEvent->InsertTrack(((AliFlowTrack*) cand));
-      }
-      if (fDebug > 0) printf("TPCevent %d\n", fFlowEvent->NumberOfTracks());
-
-      for (Int_t pTracks = 0; pTracks < unp ; pTracks++) {
-         for (Int_t nTracks = pTracks + 1; nTracks < unp ; nTracks++) {
-            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(up[pTracks], up[nTracks]))) continue;
-            if (fCandidateMinEta && (!CheckCandidateEtaPtCut(up[pTracks], up[nTracks]))) continue;
-            PtSelector(1, up[pTracks], up[nTracks]);
-         }
-      }
-      for (Int_t nTracks = 0; nTracks < unn ; nTracks++) {
-         for (Int_t pTracks = nTracks + 1; pTracks < unn ; pTracks++) {
-            if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(un[nTracks], un[pTracks]))) continue;
-            if (fCandidateEtaPtCut && (!CheckCandidateEtaPtCut(un[nTracks], un[pTracks]))) continue;
-            PtSelector(2, un[nTracks], un[pTracks]);
-         }
-      }
-      PostData(1, fOutputList);
-      PostData(2, fFlowEvent);
-   }
+}
+//_____________________________________________________________________________
+void AliAnalysisTaskPhiFlow::Exec(Option_t* c)
+{
+    // skip the event selection for SE task (e.g. for MC productions)
+    if(fSkipEventSelection) AliAnalysisTaskPhiFlow::UserExec(c);
+    // exec of task se will do event selection and call UserExec 
+    else AliAnalysisTaskSE::Exec(c);
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCandidates) const
@@ -905,7 +868,7 @@ void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCand
                 if(!fTypeMixing) { // mix with like-sign kaons
                     for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
                         if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
-                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
+//                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_up[pBuf]))) continue;
                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_up[pBuf]))) continue;
                         PtSelector(1, mix_up[pMix], buffer_up[pBuf]);
                     }
@@ -913,7 +876,7 @@ void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCand
                 else if(fTypeMixing) { // mix with unlike-sign kaons
                     for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
                         if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
-                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
+//                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_up[pMix], buffer_un[nBuf]))) continue;
                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_up[pMix], buffer_un[nBuf]))) continue;
                         PtSelector(2, mix_up[pMix], buffer_un[nBuf]);
                     }
@@ -924,7 +887,7 @@ void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCand
                 if(!fTypeMixing) { // mix with like-sign kaons
                     for(Int_t nBuf = 0; nBuf < buffer_unn; nBuf++) {
                         if(fDebug > 1 ) cout << Form(" buffered k- track %d", nBuf) << endl;
-                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
+//                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_un[nBuf]))) continue;
                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_un[nBuf]))) continue;
                         PtSelector(2, mix_un[nMix], buffer_un[nBuf]);
                     }
@@ -932,7 +895,7 @@ void AliAnalysisTaskPhiFlow::ReconstructionWithEventMixing(TObjArray* MixingCand
                 else if(fTypeMixing) { // mix with unlike-sign kaons
                     for(Int_t pBuf = 0; pBuf < buffer_unp; pBuf++) {
                         if(fDebug > 1 ) cout << Form(" buffered k+ track %d", pBuf) << endl;
-                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
+//                        if (fApplyDeltaDipCut && (!CheckDeltaDipAngle(mix_un[nMix], buffer_up[pBuf]))) continue;
                         if (fCandidateMinEta && (!CheckCandidateEtaPtCut(mix_un[nMix], buffer_up[pBuf]))) continue;
                         PtSelector(1, mix_un[nMix], buffer_up[pBuf]);
                     }
@@ -990,7 +953,8 @@ void AliAnalysisTaskPhiFlow::IsMC()
      if (fDebug>1) cout << " Received MC kaon " << endl;
      Double_t b[2] = { -99., -99.};
      Double_t bCov[3] = { -99., -99., -99.};
-     track->PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov);
+     AliAODTrack copy(*track);
+     if(!copy.PropagateToDCA(fAOD->GetPrimaryVertex(), fAOD->GetMagneticField(), 100., b, bCov)) return;
      // find corresponding mc particle
      AliAODMCParticle *partMC = (AliAODMCParticle*) arrayMC->At(TMath::Abs(track->GetLabel()));
      if (!partMC) {