]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskRhoVnModulation.cxx
from Chris
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / UserTasks / AliAnalysisTaskRhoVnModulation.cxx
index a79658dc0454d83a98944731f40f87ae6395bd3b..a2a28f87a3927c313d962397f4ce1b8b219ebc3c 100644 (file)
 #include <AliLocalRhoParameter.h>
 #include <AliAnalysisTaskRhoVnModulation.h>
 
-
 class AliAnalysisTaskRhoVnModulation;
 using namespace std;
 
 ClassImp(AliAnalysisTaskRhoVnModulation)
 
 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTaskEmcalJet("AliAnalysisTaskRhoVnModulation", kTRUE), 
-    fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(kGrid), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kTRUE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+    fDebug(0), fLocalInit(0), fAttachToEvent(kTRUE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(kGrid), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.), fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
     for(Int_t i(0); i < 10; i++) {
         fProfV2Resolution[i] = 0;
         fProfV3Resolution[i] = 0;
@@ -81,21 +80,13 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTa
         fHistRCPhiEta[i] = 0;
         fHistRhoVsRCPt[i] = 0;
         fHistRCPt[i] = 0;
-        fHistDeltaPtDeltaPhi2TPC[i] = 0;
-        fHistDeltaPtDeltaPhi2V0A[i] = 0;
-        fHistDeltaPtDeltaPhi2V0C[i] = 0;
-        fHistDeltaPtDeltaPhi3TPC[i] = 0;
-        fHistDeltaPtDeltaPhi3V0A[i] = 0;
-        fHistDeltaPtDeltaPhi3V0C[i] = 0;
+        fHistDeltaPtDeltaPhi2[i] = 0;
+        fHistDeltaPtDeltaPhi3[i] = 0;
         fHistRCPhiEtaExLJ[i] = 0;
         fHistRhoVsRCPtExLJ[i] = 0;
         fHistRCPtExLJ[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJTPC[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJV0A[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJV0C[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJTPC[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJV0A[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJV0C[i] = 0;
+        fHistDeltaPtDeltaPhi2ExLJ[i] = 0;
+        fHistDeltaPtDeltaPhi3ExLJ[i] = 0;
         /* fHistRCPhiEtaRand[i] = 0; */
         /* fHistRhoVsRCPtRand[i] = 0; */
         /* fHistRCPtRand[i] = 0; */
@@ -107,21 +98,14 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation() : AliAnalysisTa
         fHistJetPtArea[i] = 0;
         fHistJetPtConstituents[i] = 0;
         fHistJetEtaRho[i] = 0;
-        fHistJetPsiTPCPt[i] = 0;
-        fHistJetPsiVZEROAPt[i] = 0;
-        fHistJetPsiVZEROCPt[i] = 0;
-        fHistDeltaPhi2VZEROA[i] = 0;
-        fHistDeltaPhi2VZEROC[i] = 0;
-        fHistDeltaPhi2TPC[i] = 0;
-        fHistDeltaPhi3VZEROA[i] = 0;
-        fHistDeltaPhi3VZEROC[i] = 0;
-        fHistDeltaPhi3TPC[i] = 0;
+        fHistJetPsi2Pt[i] = 0;
+        fHistJetPsi3Pt[i] = 0;
    }
     // default constructor
 }
 //_____________________________________________________________________________
 AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name, runModeType type) : AliAnalysisTaskEmcalJet(name, kTRUE),
-  fDebug(0), fInitialized(0), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(1), fReduceBinsYByFactor(1), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("Q"), fRunModeType(type), fDataType(kESD), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.),  fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kTRUE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
+  fDebug(0), fLocalInit(0), fAttachToEvent(kTRUE), fFillHistograms(kTRUE), fFillQAHistograms(kTRUE), fReduceBinsXByFactor(-1.), fReduceBinsYByFactor(-1.), fNoEventWeightsForQC(kTRUE), fCentralityClasses(0), fPtBinsHybrids(0), fPtBinsJets(0), fUserSuppliedV2(0), fUserSuppliedV3(0), fUserSuppliedR2(0), fUserSuppliedR3(0), fTracksCont(0), fJetsCont(0), fUseScaledRho(0), fNAcceptedTracks(0), fNAcceptedTracksQCn(0), fFitModulationType(kNoFit), fQCRecovery(kTryFit), fUsePtWeight(kTRUE), fDetectorType(kTPC), fFitModulationOptions("QWLI"), fRunModeType(type), fDataType(kESD), fCollisionType(kPbPb), fRandom(0), fMappedRunNumber(0), fInCentralitySelection(-1), fFitModulation(0), fMinPvalue(0.01), fMaxPvalue(1), fNameJetClones(0), fNamePicoTrackClones(0), fNameRho(0), fLocalJetMinEta(-10), fLocalJetMaxEta(-10), fLocalJetMinPhi(-10), fLocalJetMaxPhi(-10), fSoftTrackMinPt(0.15), fSoftTrackMaxPt(5.),  fAbsVertexZ(10), fHistCentrality(0), fHistVertexz(0), fHistRunnumbersPhi(0), fHistRunnumbersEta(0), fHistPvaluePDF(0), fHistPvalueCDF(0), fMinDisanceRCtoLJ(0), fRandomConeRadius(-1.), fMaxCones(-1), fAbsVnHarmonics(kTRUE), fExcludeLeadingJetsFromFit(1.), fRebinSwapHistoOnTheFly(kTRUE), fPercentageOfFits(10.), fUseV0EventPlaneFromHeader(kTRUE), fExplicitOutlierCut(-1), fMinLeadingHadronPt(0), fSubtractJetPt(kFALSE), fOutputList(0), fOutputListGood(0), fOutputListBad(0), fHistAnalysisSummary(0), fHistSwap(0), fProfV2(0), fProfV2Cumulant(0), fProfV3(0), fProfV3Cumulant(0), fHistPsiControl(0), fHistPsiSpread(0), fHistPsiVZEROA(0), fHistPsiVZEROC(0), fHistPsiVZERO(0), fHistPsiTPC(0), fHistRhoVsMult(0), fHistRhoVsCent(0), fHistRhoAVsMult(0), fHistRhoAVsCent(0) {
     for(Int_t i(0); i < 10; i++) {
         fProfV2Resolution[i] = 0;
         fProfV3Resolution[i] = 0;
@@ -141,21 +125,13 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name,
         fHistRCPhiEta[i] = 0;
         fHistRhoVsRCPt[i] = 0;
         fHistRCPt[i] = 0;
-        fHistDeltaPtDeltaPhi2TPC[i] = 0;
-        fHistDeltaPtDeltaPhi2V0A[i] = 0;
-        fHistDeltaPtDeltaPhi2V0C[i] = 0;
-        fHistDeltaPtDeltaPhi3TPC[i] = 0;
-        fHistDeltaPtDeltaPhi3V0A[i] = 0;
-        fHistDeltaPtDeltaPhi3V0C[i] = 0;
+        fHistDeltaPtDeltaPhi2[i] = 0;
+        fHistDeltaPtDeltaPhi3[i] = 0;
         fHistRCPhiEtaExLJ[i] = 0;
         fHistRhoVsRCPtExLJ[i] = 0;
         fHistRCPtExLJ[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJTPC[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJV0A[i] = 0;
-        fHistDeltaPtDeltaPhi2ExLJV0C[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJTPC[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJV0A[i] = 0;
-        fHistDeltaPtDeltaPhi3ExLJV0C[i] = 0;
+        fHistDeltaPtDeltaPhi2ExLJ[i] = 0;
+        fHistDeltaPtDeltaPhi3ExLJ[i] = 0;
         /* fHistRCPhiEtaRand[i] = 0; */
         /* fHistRhoVsRCPtRand[i] = 0; */
         /* fHistRCPtRand[i] = 0; */
@@ -167,15 +143,8 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name,
         fHistJetPtArea[i] = 0;
         fHistJetPtConstituents[i] = 0;
         fHistJetEtaRho[i] = 0;
-        fHistJetPsiTPCPt[i] = 0;
-        fHistJetPsiVZEROAPt[i] = 0;
-        fHistJetPsiVZEROCPt[i] = 0;
-        fHistDeltaPhi2VZEROA[i] = 0;
-        fHistDeltaPhi2VZEROC[i] = 0;
-        fHistDeltaPhi2TPC[i] = 0;
-        fHistDeltaPhi3VZEROA[i] = 0;
-        fHistDeltaPhi3VZEROC[i] = 0;
-        fHistDeltaPhi3TPC[i] = 0;
+        fHistJetPsi2Pt[i] = 0;
+        fHistJetPsi3Pt[i] = 0;
     }
     // constructor
     DefineInput(0, TChain::Class());
@@ -188,6 +157,13 @@ AliAnalysisTaskRhoVnModulation::AliAnalysisTaskRhoVnModulation(const char* name,
         } break;
         default: fDebug = -1;   // suppress debug info explicitely when not running locally
     }
+    switch (fCollisionType) {
+        case kPythia : {
+            fFitModulationType = kNoFit;
+        } break;
+        default : break;
+    }
+    if(fLocalRhoName=="") fLocalRhoName = Form("LocalRhoFrom_%s", GetName());
 }
 //_____________________________________________________________________________
 AliAnalysisTaskRhoVnModulation::~AliAnalysisTaskRhoVnModulation()
@@ -201,14 +177,38 @@ AliAnalysisTaskRhoVnModulation::~AliAnalysisTaskRhoVnModulation()
     if(fCentralityClasses)      delete fCentralityClasses;
 }
 //_____________________________________________________________________________
+void AliAnalysisTaskRhoVnModulation::ExecOnce()
+{
+    // Init the analysis
+    fLocalRho = new AliLocalRhoParameter(fLocalRhoName.Data(), 0); 
+    if(fAttachToEvent) {
+        if(!(InputEvent()->FindListObject(fLocalRho->GetName()))) {
+            InputEvent()->AddObject(fLocalRho);
+        } else {
+            AliFatal(Form("%s: Container with name %s already present. Aborting", GetName(), fLocalRho->GetName()));
+        }
+    }
+    AliAnalysisTaskEmcalJet::ExecOnce();        // init the base class
+    AliAnalysisTaskEmcalJet::SetVzRange(-1.*fAbsVertexZ, fAbsVertexZ);
+    if(fUseScaledRho) {
+        // unscaled rho has been retrieved by the parent class, now we retrieve rho scaled
+        fRho = dynamic_cast<AliRhoParameter*>(InputEvent()->FindListObject(Form("%s_Scaled", fRho->GetName())));
+        if(!fRho) {
+            AliFatal(Form("%s: Couldn't find container for scaled rho. Aborting !", GetName()));
+        }
+    }
+    if(!GetJetContainer()) AliFatal(Form("%s: Couldn't find jet container. Aborting !", GetName()));
+}
+//_____________________________________________________________________________
 Bool_t AliAnalysisTaskRhoVnModulation::InitializeAnalysis() 
 {
     // initialize the anaysis
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(fRandomConeRadius <= 0) fRandomConeRadius = fJetRadius;
-    if(fLocalJetMinEta > -10 && fLocalJetMaxEta > -10) SetJetEtaLimits(fLocalJetMinEta, fLocalJetMaxEta);
-    if(fLocalJetMinPhi > -10 && fLocalJetMaxPhi > -10) SetJetPhiLimits(fLocalJetMinPhi, fLocalJetMaxPhi);
-    if(fMinDisanceRCtoLJ==0) fMinDisanceRCtoLJ = .5*fJetRadius;
+    if(fRandomConeRadius <= 0) fRandomConeRadius = GetJetContainer()->GetJetRadius();
+    if(fMaxCones <= 0) fMaxCones = TMath::Nint(1.8*TMath::TwoPi()/(TMath::Pi()*fRandomConeRadius*fRandomConeRadius));
+    if(fLocalJetMinEta > -10 && fLocalJetMaxEta > -10) GetJetContainer()->SetJetEtaLimits(fLocalJetMinEta, fLocalJetMaxEta);
+    if(fLocalJetMinPhi > -10 && fLocalJetMaxPhi > -10) GetJetContainer()->SetJetPhiLimits(fLocalJetMinPhi, fLocalJetMaxPhi);
+    if(fMinDisanceRCtoLJ==0) fMinDisanceRCtoLJ = .5*GetJetRadius();
     if(dynamic_cast<AliAODEvent*>(InputEvent())) fDataType = kAOD; // determine the datatype
     else if(dynamic_cast<AliESDEvent*>(InputEvent())) fDataType = kESD;
     fHistAnalysisSummary->SetBinContent(36, (int)fDataType);
@@ -243,8 +243,6 @@ Bool_t AliAnalysisTaskRhoVnModulation::InitializeAnalysis()
         case kGrid : { fFitModulationOptions += "N0"; } break;
         default : break;
     }
-    fLocalRho = new AliLocalRhoParameter(Form("local_%s", fRho->GetName()), 0); 
-    fLocalRho->SetLocalRho(fFitModulation);
     FillAnalysisSummaryHistogram();
     return kTRUE;
 }
@@ -253,7 +251,7 @@ TH1F* AliAnalysisTaskRhoVnModulation::BookTH1F(const char* name, const char* x,
 {
     // book a TH1F and connect it to the output container
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(fReduceBinsXByFactor > 0 ) bins = TMath::Nint(bins/(double)fReduceBinsXByFactor);
+    if(fReduceBinsXByFactor > 0 ) bins = TMath::Nint(bins/fReduceBinsXByFactor);
     if(!fOutputList) return 0x0;
     TString title(name);
     if(c!=-1) { // format centrality dependent histograms accordingly
@@ -271,8 +269,8 @@ TH2F* AliAnalysisTaskRhoVnModulation::BookTH2F(const char* name, const char* x,
 {
     // book a TH2F and connect it to the output container
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(fReduceBinsXByFactor > 0 ) binsx = TMath::Nint(binsx/(double)fReduceBinsXByFactor);
-    if(fReduceBinsYByFactor > 0 ) binsy = TMath::Nint(binsy/(double)fReduceBinsYByFactor);
+    if(fReduceBinsXByFactor > 0 ) binsx = TMath::Nint(binsx/fReduceBinsXByFactor);
+    if(fReduceBinsYByFactor > 0 ) binsy = TMath::Nint(binsy/fReduceBinsYByFactor);
     if(!fOutputList) return 0x0;
     TString title(name);
     if(c!=-1) { // format centrality dependent histograms accordingly
@@ -343,6 +341,7 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
     fOutputList->Add(fHistPsiSpread);
     fHistPsiVZEROA =            BookTH1F("fHistPsiVZEROA", "#Psi_{VZEROA}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
     fHistPsiVZEROC =            BookTH1F("fHistPsiVZEROC", "#Psi_{VZEROC}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
+    fHistPsiVZERO =             BookTH1F("fHistPsiVZERO", "#Psi_{VZERO}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
     fHistPsiTPC =               BookTH1F("fHistPsiTPC", "#Psi_{TPC}", 100, -.5*TMath::Pi(), .5*TMath::Pi());
     // background
     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i ++) {
@@ -354,27 +353,31 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
     fHistRhoAVsMult =           BookTH2F("fHistRhoAVsMult", "multiplicity", "#rho * A (jet) [GeV/c]", 100, 0, 4000, 100, 0, 50);
     fHistRhoAVsCent =           BookTH2F("fHistRhoAVsCent", "centrality", "#rho * A (jet) [GeV/c]", 100, 0, 100, 100, 0, 50);
 
+    TString detector("");
+    switch (fDetectorType) {
+        case kTPC : detector+="TPC";
+            break;
+        case kVZEROA : detector+="VZEROA";
+            break;
+        case kVZEROC : detector+="VZEROC";
+            break;
+        case kVZEROComb : detector+="VZEROComb";
+            break; 
+        default: break;
+    }
     // delta pt distributions
     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i ++) {
-        if(fFillQAHistograms)   fHistRCPhiEta[i] = BookTH2F("fHistRCPhiEta", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i);
+        if(fFillQAHistograms)   fHistRCPhiEta[i] = BookTH2F("fHistRCPhiEta", "#phi (RC)", "#eta (RC)", 40, 0, TMath::TwoPi(), 40, -1, 1, i);
         fHistRhoVsRCPt[i] =            BookTH2F("fHistRhoVsRCPt", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 350, i);
         fHistRCPt[i] =                 BookTH1F("fHistRCPt", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
-        if(fFillQAHistograms)   fHistRCPhiEtaExLJ[i] = BookTH2F("fHistRCPhiEtaExLJ", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i);
-        fHistDeltaPtDeltaPhi2TPC[i] =  BookTH2F("fHistDeltaPtDeltaPhi2TPC", "#phi - #Psi_{TPC}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi2V0A[i] =  BookTH2F("fHistDeltaPtDeltaPhi2V0A", "#phi - #Psi_{V0A}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi2V0C[i] =  BookTH2F("fHistDeltaPtDeltaPhi2V0C", "#phi - #Psi_{V0C}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3TPC[i] =  BookTH2F("fHistDeltaPtDeltaPhi3TPC", "#phi - #Psi_{TPC}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3V0A[i] =  BookTH2F("fHistDeltaPtDeltaPhi3V0A", "#phi - #Psi_{V0A}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3V0C[i] =  BookTH2F("fHistDeltaPtDeltaPhi3V0C", "#phi - #Psi_{V0C}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
+        if(fFillQAHistograms)   fHistRCPhiEtaExLJ[i] = BookTH2F("fHistRCPhiEtaExLJ", "#phi (RC)", "#eta (RC)", 40, 0, TMath::TwoPi(), 40, -1, 1, i);
+        fHistDeltaPtDeltaPhi2[i] =  BookTH2F("fHistDeltaPtDeltaPhi2", Form("#phi - #Psi_{2, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
+        fHistDeltaPtDeltaPhi3[i] =  BookTH2F("fHistDeltaPtDeltaPhi3", Form("#phi - #Psi_{3, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::TwoPi()/3., 400, -70, 130, i);
         fHistRhoVsRCPtExLJ[i] =        BookTH2F("fHistRhoVsRCPtExLJ", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 350, i);
         fHistRCPtExLJ[i] =             BookTH1F("fHistRCPtExLJ", "p_{t} (RC) [GeV/c]", 130, -20, 150, i);
         /* fHistRCPhiEtaRand[i] =         BookTH2F("fHistRCPhiEtaRand", "#phi (RC)", "#eta (RC)", 100, 0, TMath::TwoPi(), 100, -1, 1, i); */
-        fHistDeltaPtDeltaPhi2ExLJTPC[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJTPC", "#phi - #Psi_{TPC}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi2ExLJV0A[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJV0A", "#phi - #Psi_{V0A}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi2ExLJV0C[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJV0C", "#phi - #Psi_{V0C}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3ExLJTPC[i] = BookTH2F("fHistDeltaPtDeltaPhi3ExLJTPC", "#phi - #Psi_{TPC}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3ExLJV0A[i] = BookTH2F("fHistDeltaPtDeltaPhi3ExLJV0A", "#phi - #Psi_{V0A}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
-        fHistDeltaPtDeltaPhi3ExLJV0C[i] = BookTH2F("fHistDeltaPtDeltaPhi3ExLJV0C", "#phi - #Psi_{V0C}", "#delta p_{t} [GeV/c]", 50, 0, TMath::TwoPi()/3., 100, -50, 100, i);
+        fHistDeltaPtDeltaPhi2ExLJ[i] = BookTH2F("fHistDeltaPtDeltaPhi2ExLJ", Form("#phi - #Psi_{2, %s}", detector.Data()),  "#delta p_{t} [GeV/c]", 40, 0, TMath::Pi(), 400, -70, 130, i);
+        fHistDeltaPtDeltaPhi3ExLJ[i] = BookTH2F("fHistDeltaPtDeltaPhi3ExLJ", Form("#phi - #Psi_{3, %s}", detector.Data()), "#delta p_{t} [GeV/c]", 40, 0, TMath::TwoPi()/3., 400, -70, 130, i);
         /* fHistRhoVsRCPtRand[i] =        BookTH2F("fHistRhoVsRCPtRand", "p_{t} (RC) [GeV/c]", "#rho * A (RC) [GeV/c]", 100, 0, 300, 100, 0, 350, i); */
         /* fHistRCPtRand[i] =             BookTH1F("fHistRCPtRand", "p_{t} (RC) [GeV/c]", 130, -20, 150, i); */
         /* fHistDeltaPtDeltaPhi2Rand[i] =  BookTH2F("fHistDeltaPtDeltaPhi2Rand", "#phi - #Psi_{TPC}", "#delta p_{t} [GeV/c]", 50, 0, TMath::Pi(), 100, -50, 100, i); */
@@ -387,32 +390,30 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
         fHistJetPtConstituents[i] =    BookTH2F("fHistJetPtConstituents", "p_{t} [GeV/c]", "Area", 350, -100, 250, 60, 0, 150, i);
         fHistJetEtaRho[i] =            BookTH2F("fHistJetEtaRho", "#eta", "#rho", 100, -1, 1, 100, 0, 300, i);
         // in plane and out of plane spectra
-        fHistJetPsiTPCPt[i] =          BookTH2F("fHistJetPsiTPCPt", "#phi_{jet} - #Psi_{2, TPC}", "p_{t} [GeV/c]", 50, 0., TMath::Pi(), 700, -100, 250, i);
-        fHistJetPsiVZEROAPt[i] =       BookTH2F("fHistJetPsiVZEROAPt", "#phi_{jet} - #Psi_{2, VZEROA}", "p_{t} [GeV/c]", 50, 0., TMath::Pi(), 700, -100, 250, i);
-        fHistJetPsiVZEROCPt[i] =       BookTH2F("fHistJetPsiVZEROCPt", "#phi_{jet} - #Psi_{V2, ZEROC}", "p_{t} [GeV/c]", 50, 0., TMath::Pi(), 700, -100, 250, i);
-        // phi minus psi
-        fHistDeltaPhi2VZEROA[i] =       BookTH1F("fHistDeltaPhi2VZEROA", "#phi_{jet} - #Psi_{2, VZEROA}", 50, 0, TMath::Pi(), i);
-        fHistDeltaPhi2VZEROC[i] =       BookTH1F("fHistDeltaPhi2VZEROC", "#phi_{jet} - #Psi_{2, VZEROC}", 50, 0, TMath::Pi(), i);
-        fHistDeltaPhi2TPC[i] =          BookTH1F("fHistDeltaPhi2TPC", "#phi_{jet} - #Psi_{2, TPC}", 50, 0, TMath::Pi(), i);
-        fHistDeltaPhi3VZEROA[i] =       BookTH1F("fHistDeltaPhi3VZEROA", "#phi_{jet} - #Psi_{2, VZEROA}", 50, 0, TMath::TwoPi()/3., i);
-        fHistDeltaPhi3VZEROC[i] =       BookTH1F("fHistDeltaPhi3VZEROC", "#phi_{jet} - #Psi_{2, VZEROC}", 50, 0, TMath::TwoPi()/3., i);
-        fHistDeltaPhi3TPC[i] =          BookTH1F("fHistDeltaPhi3TPC", "#phi_{jet} - #Psi_{2, TPC}", 50, 0, TMath::TwoPi()/3., i);
-
-        fProfV2Resolution[i] = new TProfile(Form("fProfV2Resolution_%i", i), Form("fProfV2Resolution_%i", i), 8, -0.5, 7.5);
+        fHistJetPsi2Pt[i] =          BookTH2F("fHistJetPsi2Pt", Form("#phi_{jet} - #Psi_{2, %s}", detector.Data()), "p_{t} [GeV/c]", 40, 0., TMath::Pi(), 350, -100, 250, i);
+        fHistJetPsi3Pt[i] =          BookTH2F("fHistJetPsi3Pt", Form("#phi_{jet} - #Psi_{3, %s}", detector.Data()), "p_{t} [GeV/c]", 40, 0., TMath::TwoPi()/3., 350, -100, 250, i);
+        // profiles for all correlator permutations which are necessary to calculate each second and third order event plane resolution
+        fProfV2Resolution[i] = new TProfile(Form("fProfV2Resolution_%i", i), Form("fProfV2Resolution_%i", i), 11, -0.5, 10.5);
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(2(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(2(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(2(#Psi_{VZEROA} - #Psi_{TPC}))>");
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(2(#Psi_{TPC} - #Psi_{VZEROA}))>");
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(2(#Psi_{VZEROC} - #Psi_{TPC}))>");
         fProfV2Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(2(#Psi_{TPC} - #Psi_{VZEROC}))>");
+        fProfV2Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_A}))>");
+        fProfV2Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(2(#Psi_{VZERO} - #Psi_{TPC_B}))>");
+        fProfV2Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(2(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
         fOutputList->Add(fProfV2Resolution[i]); 
-        fProfV3Resolution[i] = new TProfile(Form("fProfV3Resolution_%i", i), Form("fProfV3Resolution_%i", i), 8, -0.5, 7.5);
+        fProfV3Resolution[i] = new TProfile(Form("fProfV3Resolution_%i", i), Form("fProfV3Resolution_%i", i), 11, -0.5, 10.5);
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(3, "<cos(3(#Psi_{VZEROA} - #Psi_{VZEROC}))>");
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(4, "<cos(3(#Psi_{VZEROC} - #Psi_{VZEROA}))>");
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(5, "<cos(3(#Psi_{VZEROA} - #Psi_{TPC}))>");
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(6, "<cos(3(#Psi_{TPC} - #Psi_{VZEROA}))>");
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(7, "<cos(3(#Psi_{VZEROC} - #Psi_{TPC}))>");
         fProfV3Resolution[i]->GetXaxis()->SetBinLabel(8, "<cos(3(#Psi_{TPC} - #Psi_{VZEROC}))>");
+        fProfV3Resolution[i]->GetXaxis()->SetBinLabel(9, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_A}))>");
+        fProfV3Resolution[i]->GetXaxis()->SetBinLabel(10, "<cos(3(#Psi_{VZERO} - #Psi_{TPC_B}))>");
+        fProfV3Resolution[i]->GetXaxis()->SetBinLabel(11, "<cos(3(#Psi_{TPC_A} - #Psi_{TPC_B}))>");
         fOutputList->Add(fProfV3Resolution[i]); 
     }
     // cdf and pdf of chisquare distribution
@@ -451,7 +452,7 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
         fHistRunnumbersPhi->Sumw2();
         fOutputList->Add(fHistRunnumbersPhi);
     }
-    fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 50, -0.5, 50.5);
+    fHistAnalysisSummary = BookTH1F("fHistAnalysisSummary", "flag", 52, -0.5, 52.5);
     fHistSwap = new TH1F("fHistSwap", "fHistSwap", 20, 0, TMath::TwoPi());
     if(fUsePtWeight) fHistSwap->Sumw2();
 
@@ -474,14 +475,18 @@ void AliAnalysisTaskRhoVnModulation::UserCreateOutputObjects()
         } break;
         default: break;
     }
+
+    // get the containers
+    fTracksCont = GetParticleContainer("Tracks");
+    fJetsCont = GetJetContainer("Jets");
 }
 //_____________________________________________________________________________
 Bool_t AliAnalysisTaskRhoVnModulation::Run()
 {
     // user exec: execute once for each event
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(!(fTracks||fJets||fRho)) return kFALSE;
-    if(!fInitialized) fInitialized = InitializeAnalysis();
+    if(!fTracks||!fJets||!fRho) return kFALSE;
+    if(!fLocalInit) fLocalInit = InitializeAnalysis();
     // reject the event if expected data is missing
     if(!PassesCuts(InputEvent())) return kFALSE;
     if(!fCaloClusters && fDebug > 0) printf(" > Warning: couldn't retreive calo clusters! < \n");
@@ -509,7 +514,17 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
         default : break;
     }
     switch (fFitModulationType) { // do the fits
-        case kNoFit : { fFitModulation->FixParameter(0, fLocalRho->GetVal()); } break;
+        case kNoFit : { 
+             switch (fCollisionType) {
+                 case kPythia : { // background is zero for pp jets
+                     fFitModulation->FixParameter(0, 0);
+                     fLocalRho->SetVal(0);
+                 } break;
+                 default :  {
+                     fFitModulation->FixParameter(0, fLocalRho->GetVal()); 
+                 } break;
+             }
+        } break;
         case kV2 : {    // only v2
             if(CorrectRho(psi2, psi3)) {
                 fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
@@ -517,7 +532,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     Double_t r(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
                     if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
                 }
-                CalculateEventPlaneResolution(vzero, tpc);
+                CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
             }
         } break;
         case kV3 : {    // only v3
@@ -527,7 +542,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     if(r > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r);
                 }
                 fProfV3->Fill(fCent, fFitModulation->GetParameter(3));
-                CalculateEventPlaneResolution(vzero, tpc);
+                CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
             }
         } break;
         case kQC2 : {   // qc2 analysis
@@ -537,7 +552,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
                     Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
                     if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
-                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
+                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
                 }
                 if (fUsePtWeight) { // use weighted weights
                     Double_t dQCnM11 = (fNoEventWeightsForQC) ? 1. : QCnM11();
@@ -548,7 +563,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     fProfV2->Fill(fCent, fFitModulation->GetParameter(3), dQCnM*(dQCnM-1));
                     fProfV3->Fill(fCent, fFitModulation->GetParameter(7), dQCnM*(dQCnM-1));
                 }
-                CalculateEventPlaneResolution(vzero, tpc);
+                CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
             }
         } break;
         case kQC4 : {
@@ -558,7 +573,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
                     Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
                     if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)*r2);
-                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)*r3);
+                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)*r3);
                 }
                 if (fUsePtWeight) { // use weighted weights
                     fProfV2->Fill(fCent, TMath::Power(fFitModulation->GetParameter(3),0.5)/*, QCnM1111()*/);
@@ -568,7 +583,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     fProfV3->Fill(fCent, TMath::Power(fFitModulation->GetParameter(7),0.5)/*, QCnM()*(QCnM()-1)*(QCnM()-2)*(QCnM()-3)*/);
                 }
             }
-            CalculateEventPlaneResolution(vzero, tpc);
+            CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
         } break;
         default : {
             if(CorrectRho(psi2, psi3)) {
@@ -576,16 +591,19 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
                     Double_t r2(fUserSuppliedR2->GetBinContent(fUserSuppliedR2->GetXaxis()->FindBin(fCent)));
                     Double_t r3(fUserSuppliedR3->GetBinContent(fUserSuppliedR3->GetXaxis()->FindBin(fCent)));
                     if(r2 > 0) fFitModulation->SetParameter(3, fFitModulation->GetParameter(3)/r2);
-                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(3)/r3);
+                    if(r3 > 0) fFitModulation->SetParameter(7, fFitModulation->GetParameter(7)/r3);
                 }
                 fProfV2->Fill(fCent, fFitModulation->GetParameter(3));
                 fProfV3->Fill(fCent, fFitModulation->GetParameter(7));
-                CalculateEventPlaneResolution(vzero, tpc);
+                CalculateEventPlaneResolution(vzero, vzeroComb, tpc);
             }
         } break;
     }
+    // if all went well, update the local rho parameter
+    fLocalRho->SetLocalRho(fFitModulation);
     // fill a number of histograms 
-    FillHistogramsAfterSubtraction(vzero, tpc);
+    if(fFillHistograms)         FillHistogramsAfterSubtraction(psi2, psi3, vzero, vzeroComb, tpc);
+    if(fFillQAHistograms)       FillQAHistograms(InputEvent());
     // send the output to the connected output container
     PostData(1, fOutputList);
     switch (fRunModeType) {
@@ -595,6 +613,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::Run()
         } break;
         default: break;
     }
+
     return kTRUE;
 }
 //_____________________________________________________________________________
@@ -643,26 +662,16 @@ void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneTPC(Double_t* tpc)
    Double_t qx2(0), qy2(0);     // for psi2
    Double_t qx3(0), qy3(0);     // for psi3
    if(fTracks) {
-       Float_t excludeInEta[] = {-999, -999};
+       Float_t excludeInEta = -999;
        if(fExcludeLeadingJetsFromFit > 0 ) {    // remove the leading jet from ep estimate
-           AliEmcalJet* leadingJet[] = {0x0, 0x0};
-           static Int_t lJets[9999] = {-1};
-           GetSortedArray(lJets, fJets);
-           for(Int_t i(0); i < fJets->GetEntriesFast(); i++) {     // get the two leading jets
-               if (1 + i > fJets->GetEntriesFast()) break;
-               leadingJet[0] = static_cast<AliEmcalJet*>(fJets->At(lJets[i]));
-               leadingJet[1] = static_cast<AliEmcalJet*>(fJets->At(lJets[i+1]));
-               if(PassesCuts(leadingJet[0]) && PassesCuts(leadingJet[1])) break;
-           }
-           if(leadingJet[0] && leadingJet[1]) {
-               for(Int_t i(0); i < 2; i++) excludeInEta[i] = leadingJet[i]->Eta();
-           }
+           AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
+           if(leadingJet) excludeInEta = leadingJet->Eta();
        }
        Int_t iTracks(fTracks->GetEntriesFast());
        for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
            AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
            if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
-           if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta[0]) < fJetRadius*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - fJetRadius - fJetMaxEta ) > 0 )) continue;
+           if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
            fNAcceptedTracks++;
            qx2+= TMath::Cos(2.*track->Phi());
            qy2+= TMath::Sin(2.*track->Phi());
@@ -677,16 +686,16 @@ void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneTPC(Double_t* tpc)
 void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneCombinedVZERO(Double_t* comb) const
 {
     // grab the combined vzero event plane
-    if(fUseV0EventPlaneFromHeader) {    // use the vzero from the header
+//    if(fUseV0EventPlaneFromHeader) {    // use the vzero from the header
         Double_t a(0), b(0), c(0), d(0);
         comb[0] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, a, b);
-        comb[1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 2, c, d);
-    } else {
-        Double_t qx2a(0), qy2a(0), qx2c(0), qy2c(0), qx3a(0), qy3a(0), qx3c(0), qy3c(0);
-        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, qx2a, qy2a);
-        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, qx2c, qy2c);
-        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, qx3a, qy3a);
-        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, qx3c, qy3c);
+        comb[1] = InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 10, 3, c, d);
+//    } else {
+//        Double_t qx2a(0), qy2a(0), qx2c(0), qy2c(0), qx3a(0), qy3a(0), qx3c(0), qy3c(0);
+//        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 2, qx2a, qy2a);
+//        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 2, qx2c, qy2c);
+//        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 8, 3, qx3a, qy3a);
+//        InputEvent()->GetEventplane()->CalculateVZEROEventPlane(InputEvent(), 9, 3, qx3c, qy3c);
 // FIXME the rest of this function isn't impelmented yet (as of 01-07-2013)
 //        Double_t chi2A(-1), chi2C(-1), chi3A(-1), chi3C(-1);     // get chi from the resolution
 //        Double_t qx2(chi2A*chi2A*qx2a+chi2C*chi2C*qx2c);
@@ -695,10 +704,10 @@ void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneCombinedVZERO(Double_t*
 //        Double_t qy3(chi3A*chi3A*qy3a+chi3C*chi3C*qy3c);
 //        comb[0] = .5*TMath::ATan2(qy2, qx2);
 //        comb[1] = (1./3.)*TMath::ATan2(qy3, qx3);
-    }
+//    }
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneResolution(Double_t vzero[2][2], Double_t* tpc) const
+void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneResolution(Double_t vzero[2][2], Double_t* vzeroComb, Double_t* tpc)
 {
     // fill the profiles for the resolution parameters
     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
@@ -714,8 +723,40 @@ void AliAnalysisTaskRhoVnModulation::CalculateEventPlaneResolution(Double_t vzer
     fProfV3Resolution[fInCentralitySelection]->Fill(5., TMath::Cos(3.*(tpc[0] - vzero[0][0])));
     fProfV3Resolution[fInCentralitySelection]->Fill(6., TMath::Cos(3.*(vzero[1][0] - tpc[0])));
     fProfV3Resolution[fInCentralitySelection]->Fill(7., TMath::Cos(3.*(tpc[0] - vzero[1][0])));
-    // FIXME no VZEROComb resolution available (as of 01-07-2013)
-}
+    // for the resolution of the combined vzero event plane, use two tpc halves as uncorrelated subdetectors
+    Double_t qx2a(0), qy2a(0);     // for psi2a, negative eta
+    Double_t qx3a(0), qy3a(0);     // for psi3a, negative eta
+    Double_t qx2b(0), qy2b(0);     // for psi2a, positive eta
+    Double_t qx3b(0), qy3b(0);     // for psi3a, positive eta
+    if(fTracks) {
+       Int_t iTracks(fTracks->GetEntriesFast());
+       for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
+           AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
+           if(!PassesCuts(track) || track->Pt() < fSoftTrackMinPt || track->Pt() > fSoftTrackMaxPt) continue;
+           if(track->Eta() < 0 ) {
+               qx2a+= TMath::Cos(2.*track->Phi());
+               qy2a+= TMath::Sin(2.*track->Phi());
+               qx3a+= TMath::Cos(3.*track->Phi());
+               qy3a+= TMath::Sin(3.*track->Phi());
+           } else if (track->Eta() > 0) {
+               qx2b+= TMath::Cos(2.*track->Phi());
+               qy2b+= TMath::Sin(2.*track->Phi());
+               qx3b+= TMath::Cos(3.*track->Phi());
+               qy3b+= TMath::Sin(3.*track->Phi());
+           }
+       }
+   }
+   Double_t tpca2(.5*TMath::ATan2(qy2a, qx2a));
+   Double_t tpca3((1./3.)*TMath::ATan2(qy3a, qx3a));
+   Double_t tpcb2(.5*TMath::ATan2(qy2b, qx2b));
+   Double_t tpcb3((1./3.)*TMath::ATan2(qy3b, qx3b));
+   fProfV2Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(2.*(vzeroComb[0] - tpca2)));
+   fProfV2Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(2.*(vzeroComb[0] - tpcb2)));
+   fProfV2Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(2.*(tpca2 - tpcb2))); 
+   fProfV3Resolution[fInCentralitySelection]->Fill(8., TMath::Cos(3.*(vzeroComb[1] - tpca3)));
+   fProfV3Resolution[fInCentralitySelection]->Fill(9., TMath::Cos(3.*(vzeroComb[1] - tpcb3)));
+   fProfV3Resolution[fInCentralitySelection]->Fill(10., TMath::Cos(3.*(tpca3 - tpcb3))); 
+}   
 //_____________________________________________________________________________
 Double_t AliAnalysisTaskRhoVnModulation::CalculateEventPlaneChi(Double_t resEP) const
 {
@@ -729,7 +770,7 @@ Double_t AliAnalysisTaskRhoVnModulation::CalculateEventPlaneChi(Double_t resEP)
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskRhoVnModulation::CalculateRandomCone(Float_t &pt, Float_t &eta, Float_t &phi, 
-        AliEmcalJet* jet, Bool_t randomize) const
+        AliEmcalJet* jet) const
 {
     // get a random cone
     if(fDebug > 1) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
@@ -740,15 +781,15 @@ void AliAnalysisTaskRhoVnModulation::CalculateRandomCone(Float_t &pt, Float_t &e
         phiJet = jet->Phi();
     }
     // force the random cones to at least be within detector acceptance
-    Float_t minPhi(fJetMinPhi), maxPhi(fJetMaxPhi);
+    Float_t minPhi(GetJetContainer()->GetJetPhiMin()), maxPhi(GetJetContainer()->GetJetPhiMax());
     if(maxPhi > TMath::TwoPi()) maxPhi = TMath::TwoPi();
     if(minPhi < 0 ) minPhi = 0;
-    Float_t diffRcRJR(TMath::Abs(fRandomConeRadius-fJetRadius));
+    Float_t diffRcRJR(TMath::Abs(fRandomConeRadius-GetJetContainer()->GetJetRadius()));
     // construct a random cone and see if it's far away enough from the leading jet
     Int_t attempts(1000);
     while(kTRUE) {
         attempts--;
-        eta = gRandom->Uniform(fJetMinEta+diffRcRJR, fJetMaxEta-diffRcRJR);
+        eta = gRandom->Uniform(GetJetContainer()->GetJetEtaMin()+diffRcRJR, GetJetContainer()->GetJetEtaMax()-diffRcRJR);
         phi = gRandom->Uniform(minPhi, maxPhi);
 
         dJet = TMath::Sqrt((etaJet-eta)*(etaJet-eta)+(phiJet-phi)*(phiJet-phi));
@@ -758,21 +799,15 @@ void AliAnalysisTaskRhoVnModulation::CalculateRandomCone(Float_t &pt, Float_t &e
             return;
         }
     }
-    if(fTracks) {
-        Int_t iTracks(fTracks->GetEntriesFast());
-        for(Int_t i(0); i < iTracks; i++) {
-            AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
-            if(!PassesCuts(track)) continue;
+    if(fTracksCont) {
+        AliVParticle* track = fTracksCont->GetNextAcceptParticle(0);
+        while(track) {
             Float_t etaTrack(track->Eta()), phiTrack(track->Phi()), ptTrack(track->Pt());
-            // if requested, randomize eta and phi to destroy any correlated fluctuations
-            if(randomize) {
-                etaTrack = gRandom->Uniform(fTrackMinEta, fTrackMaxEta);
-                phiTrack = gRandom->Uniform(minPhi, maxPhi);
-            }
             // get distance from cone
             if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi + TMath::TwoPi())) phiTrack+=TMath::TwoPi();
             if(TMath::Abs(phiTrack-phi) > TMath::Abs(phiTrack - phi - TMath::TwoPi())) phiTrack-=TMath::TwoPi();
             if(TMath::Sqrt(TMath::Abs((etaTrack-eta)*(etaTrack-eta)+(phiTrack-phi)*(phiTrack-phi))) <= fRandomConeRadius) pt+=ptTrack;
+            track = fTracksCont->GetNextAcceptParticle();
         }
     }
 }
@@ -871,7 +906,7 @@ void AliAnalysisTaskRhoVnModulation::QCnDiffentialFlowVectors(
             for(Int_t ptBin(0); ptBin < ptBins->GetSize()-1; ptBin++) {
                 AliEmcalJet* poi = static_cast<AliEmcalJet*>(pois->At(i));
                 if(PassesCuts(poi)) {    
-                    Double_t pt(poi->Pt()-poi->Area()*fLocalRho->GetLocalVal(poi->Phi(), fJetRadius, fLocalRho->GetVal()));
+                    Double_t pt(poi->Pt()-poi->Area()*fLocalRho->GetLocalVal(poi->Phi(), GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
                     if(pt >= ptBins->At(ptBin) && pt < ptBins->At(ptBin+1)) {    
                             repn[ptBin]+=TMath::Cos(((double)n)*poi->Phi());
                             impn[ptBin]+=TMath::Sin(((double)n)*poi->Phi());
@@ -1054,26 +1089,16 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
         default: break;
     }
     Int_t iTracks(fTracks->GetEntriesFast());
-    Double_t excludeInEta[] = {-999, -999};
-    Double_t excludeInPhi[] = {-999, -999};
-    Double_t excludeInPt[]  = {-999, -999};
+    Double_t excludeInEta = -999;
+    Double_t excludeInPhi = -999;
+    Double_t excludeInPt  = -999;
     if(iTracks <= 0 || fLocalRho->GetVal() <= 0 ) return kFALSE;   // no use fitting an empty event ...
     if(fExcludeLeadingJetsFromFit > 0 ) {
-        AliEmcalJet* leadingJet[] = {0x0, 0x0};
-        static Int_t lJets[9999] = {-1};
-        GetSortedArray(lJets, fJets);
-        for(Int_t i(0); i < fJets->GetEntriesFast(); i++) {     // get the two leading jets
-            if (1 + i > fJets->GetEntriesFast()) break;
-            leadingJet[0] = static_cast<AliEmcalJet*>(fJets->At(lJets[i]));
-            leadingJet[1] = static_cast<AliEmcalJet*>(fJets->At(lJets[i+1]));
-            if(PassesCuts(leadingJet[0]) && PassesCuts(leadingJet[1])) break;
-        }
-        if(leadingJet[0] && leadingJet[1]) {
-            for(Int_t i(0); i < 2; i++) {
-                excludeInEta[i] = leadingJet[i]->Eta();
-                excludeInPhi[i] = leadingJet[i]->Phi();
-                excludeInPt[i]  = leadingJet[i]->Pt();
-            }
+        AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
+        if(leadingJet) {
+            excludeInEta = leadingJet->Eta();
+            excludeInPhi = leadingJet->Phi();
+            excludeInPt = leadingJet->Pt();
         }
     }
     fHistSwap->Reset();                 // clear the histogram
@@ -1081,11 +1106,12 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
     if(fRebinSwapHistoOnTheFly) {
         if(fNAcceptedTracks < 49) fNAcceptedTracks = 49;       // avoid aliasing effects
         _tempSwap = TH1F("_tempSwap", "_tempSwap", TMath::CeilNint(TMath::Sqrt(fNAcceptedTracks)), 0, TMath::TwoPi());
+        if(fUsePtWeight) _tempSwap.Sumw2();
     }
     else _tempSwap = *fHistSwap;         // now _tempSwap holds the desired histo
     for(Int_t i(0); i < iTracks; i++) {
             AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(i));
-            if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta[0]) < fJetRadius*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - fJetRadius - fJetMaxEta ) > 0 )) continue;
+            if(fExcludeLeadingJetsFromFit > 0 &&( (TMath::Abs(track->Eta() - excludeInEta) < GetJetContainer()->GetJetRadius()*fExcludeLeadingJetsFromFit ) || (TMath::Abs(track->Eta()) - GetJetContainer()->GetJetRadius() - GetJetContainer()->GetJetEtaMax() ) > 0 )) continue;
             if(!PassesCuts(track) || track->Pt() > fSoftTrackMaxPt || track->Pt() < fSoftTrackMinPt) continue;
             if(fUsePtWeight) _tempSwap.Fill(track->Phi(), track->Pt());
             else _tempSwap.Fill(track->Phi());
@@ -1152,12 +1178,8 @@ Bool_t AliAnalysisTaskRhoVnModulation::CorrectRho(Double_t psi2, Double_t psi3)
                 }
                 if(fExcludeLeadingJetsFromFit) {       // visualize the excluded region
                     TF2 *f2 = new TF2(Form("%s_LJ", didacticSurface->GetName()),"[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", 0, TMath::TwoPi(), -1, 1);
-                    f2->SetParameters(excludeInPt[0]/3.,excludeInPhi[0],.1,excludeInEta[0],.1);
+                    f2->SetParameters(excludeInPt/3.,excludeInPhi,.1,excludeInEta,.1);
                     didacticSurface->GetListOfFunctions()->Add(f2);
-                    TF2 *f3 = new TF2(Form("%s_NLJ", didacticSurface->GetName()),"[0]*TMath::Gaus(x,[1],[2])*TMath::Gaus(y,[3],[4])", 0, TMath::TwoPi(), -1, 1);
-                    f3->SetParameters(excludeInPt[1]/3.,excludeInPhi[1],.1,excludeInEta[1],.1);
-                    f3->SetLineColor(kGreen);
-                    didacticSurface->GetListOfFunctions()->Add(f3);
                 }
                 fOutputListGood->Add(didacticSurface);
             } break;
@@ -1194,7 +1216,7 @@ Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(AliVEvent* event)
 {
     // event cuts
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    if(!event) return kFALSE;
+    if(!event || !AliAnalysisTaskEmcal::IsEventSelected()) return kFALSE;
     if(TMath::Abs(InputEvent()->GetPrimaryVertex()->GetZ()) > 10.) return kFALSE;
     // aod and esd specific checks
     switch (fDataType) {
@@ -1211,15 +1233,18 @@ Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(AliVEvent* event)
     fCent = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
     if(fCent <= fCentralityClasses->At(0) || fCent >= fCentralityClasses->At(fCentralityClasses->GetSize()-1) || TMath::Abs(fCent-InputEvent()->GetCentrality()->GetCentralityPercentile("TRK")) > 5.) return kFALSE;
     // determine centrality class
+    fInCentralitySelection = -1;
     for(Int_t i(0); i < fCentralityClasses->GetSize()-1; i++) {
         if(fCent >= fCentralityClasses->At(i) && fCent <= fCentralityClasses->At(1+i)) {
             fInCentralitySelection = i;
             break; }
     } 
+    if(fInCentralitySelection<0) return kFALSE;     // should be null op
     if(fExplicitOutlierCut == 2010 || fExplicitOutlierCut == 2011) {
        if(!PassesCuts(fExplicitOutlierCut)) return kFALSE;
     }
-    if(fFillQAHistograms) FillQAHistograms(event);
+    if(fRho->GetVal() <= 0 ) return kFALSE;
+    if(fTracks->GetEntries() < 1) return kFALSE;
     return kTRUE;
 }
 //_____________________________________________________________________________
@@ -1253,18 +1278,17 @@ Bool_t AliAnalysisTaskRhoVnModulation::PassesCuts(const AliVCluster* cluster) co
     return kTRUE;
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillHistogramsAfterSubtraction(Double_t vzero[2][2], Double_t* tpc) const
+void AliAnalysisTaskRhoVnModulation::FillHistogramsAfterSubtraction(Double_t psi2, Double_t psi3, Double_t vzero[2][2], Double_t* vzeroComb, Double_t* tpc)
 {
     // fill histograms 
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
     FillTrackHistograms();
     /* FillClusterHistograms(); */
-    FillJetHistograms(vzero, tpc); 
+    FillJetHistograms(psi2, psi3); 
     /* FillCorrectedClusterHistograms(); */
-    FillEventPlaneHistograms(vzero, tpc);
+    FillEventPlaneHistograms(vzero, vzeroComb, tpc);
     FillRhoHistograms();
-    FillDeltaPtHistograms(vzero, tpc);
-    FillDeltaPhiHistograms(vzero, tpc);
+    FillDeltaPtHistograms(psi2, psi3);
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskRhoVnModulation::FillTrackHistograms() const
@@ -1305,7 +1329,7 @@ void AliAnalysisTaskRhoVnModulation::FillCorrectedClusterHistograms() const
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillEventPlaneHistograms(Double_t vzero[2][2], Double_t* tpc) const
+void AliAnalysisTaskRhoVnModulation::FillEventPlaneHistograms(Double_t vzero[2][2], Double_t* vzeroComb, Double_t* tpc) const
 {
     // fill event plane histograms
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
@@ -1317,89 +1341,65 @@ void AliAnalysisTaskRhoVnModulation::FillEventPlaneHistograms(Double_t vzero[2][
     fHistPsiControl->Fill(7.5, tpc[1]);         // tpc psi 3
     fHistPsiVZEROA->Fill(vzero[0][0]);
     fHistPsiVZEROC->Fill(vzero[1][0]);
+    fHistPsiVZERO->Fill(vzeroComb[0]);
     fHistPsiTPC->Fill(tpc[0]);
     fHistPsiSpread->Fill(0.5, TMath::Abs(vzero[0][0]-vzero[1][0]));
     fHistPsiSpread->Fill(1.5, TMath::Abs(vzero[0][0]-tpc[0]));
     fHistPsiSpread->Fill(2.5, TMath::Abs(vzero[1][0]-tpc[0]));
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillRhoHistograms() const
+void AliAnalysisTaskRhoVnModulation::FillRhoHistograms()
 {
     // fill rho histograms
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
     fHistRhoPackage[fInCentralitySelection]->Fill(fLocalRho->GetVal());    // save the rho estimate from the emcal jet package
     // get multiplicity FIXME inefficient
-    Int_t iTracks(fTracks->GetEntriesFast()), mult(0), iJets(fJets->GetEntriesFast());
-    for(Int_t i(0); i < iTracks; i ++) { if(PassesCuts(static_cast<AliVTrack*>(fTracks->At(i)))) mult++; }
+    Int_t iJets(fJets->GetEntriesFast());
     Double_t rho(fLocalRho->GetLocalVal(TMath::Pi(), TMath::Pi(), fLocalRho->GetVal()));
     fHistRho[fInCentralitySelection]->Fill(rho);
-    fHistRhoVsMult->Fill(mult, rho);
+    fHistRhoVsMult->Fill(fTracks->GetEntries(), rho);
     fHistRhoVsCent->Fill(fCent, rho);
     for(Int_t i(0); i < iJets; i++) {
         AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
         if(!PassesCuts(jet)) continue;
-        fHistRhoAVsMult->Fill(mult, rho * jet->Area());
+        fHistRhoAVsMult->Fill(fTracks->GetEntries(), rho * jet->Area());
         fHistRhoAVsCent->Fill(fCent, rho * jet->Area());
     }
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillDeltaPtHistograms(Double_t vzero[2][2], Double_t* tpc) const
+void AliAnalysisTaskRhoVnModulation::FillDeltaPtHistograms(Double_t psi2, Double_t psi3) const
 {
     // fill delta pt histograms
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    Int_t i(0), maxCones(20);
-    AliEmcalJet* leadingJet(0x0);
-    static Int_t sJets[9999] = {-1};
-    GetSortedArray(sJets, fJets);
-    do { // get the leading jet 
-        leadingJet = static_cast<AliEmcalJet*>(fJets->At(sJets[i]));
-        i++;
-    }
-    while (!PassesCuts(leadingJet)&&i<fJets->GetEntriesFast()); 
+    Int_t i(0);
+    AliEmcalJet* leadingJet(GetJetContainer()->GetLeadingJet());
     if(!leadingJet && fDebug > 0) printf(" > failed to retrieve leading jet ! < \n");
     const Float_t areaRC = fRandomConeRadius*fRandomConeRadius*TMath::Pi();
     // we're retrieved the leading jet, now get a random cone
-    for(i = 0; i < maxCones; i++) {
+    for(i = 0; i < fMaxCones; i++) {
        Float_t pt(0), eta(0), phi(0);
        // get a random cone without constraints on leading jet position
        CalculateRandomCone(pt, eta, phi, 0x0);
        if(pt > 0) {
            if(fFillQAHistograms) fHistRCPhiEta[fInCentralitySelection]->Fill(phi, eta);
-           fHistRhoVsRCPt[fInCentralitySelection]->Fill(pt, fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal())*areaRC);
+           fHistRhoVsRCPt[fInCentralitySelection]->Fill(pt, fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal())*areaRC);
            fHistRCPt[fInCentralitySelection]->Fill(pt);
-           fHistDeltaPtDeltaPhi2TPC[fInCentralitySelection]->Fill(PhaseShift(phi-tpc[0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi2V0A[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[0][0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi2V0C[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[1][0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3TPC[fInCentralitySelection]->Fill(PhaseShift(phi-tpc[1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3V0A[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[0][1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3V0C[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[1][1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
+           fHistDeltaPtDeltaPhi2[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
+           fHistDeltaPtDeltaPhi3[fInCentralitySelection]->Fill(PhaseShift(phi-psi3, 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
        }
        // get a random cone excluding leading jet area
        CalculateRandomCone(pt, eta, phi, leadingJet);
        if(pt > 0) {
            if(fFillQAHistograms) fHistRCPhiEtaExLJ[fInCentralitySelection]->Fill(phi, eta);
-           fHistRhoVsRCPtExLJ[fInCentralitySelection]->Fill(pt, fLocalRho->GetLocalVal(phi, fJetRadius, fRho->GetVal())*areaRC);
+           fHistRhoVsRCPtExLJ[fInCentralitySelection]->Fill(pt, fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal())*areaRC);
            fHistRCPtExLJ[fInCentralitySelection]->Fill(pt);
-           fHistDeltaPtDeltaPhi2ExLJTPC[fInCentralitySelection]->Fill(PhaseShift(phi-tpc[0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi2ExLJV0A[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[0][0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi2ExLJV0C[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[1][0], 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3ExLJTPC[fInCentralitySelection]->Fill(PhaseShift(phi-tpc[1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3ExLJV0A[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[0][1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
-           fHistDeltaPtDeltaPhi3ExLJV0C[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[1][1], 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
+           fHistDeltaPtDeltaPhi2ExLJ[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt - areaRC*fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
+           fHistDeltaPtDeltaPhi3ExLJ[fInCentralitySelection]->Fill(PhaseShift(phi-psi3, 3.), pt - areaRC*fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
        }
-       // get a random cone in an event with randomized phi and eta
-       /* CalculateRandomCone(pt, eta, phi, 0x0, kTRUE);
-       if( pt > 0) {
-           fHistRCPhiEtaRand[fInCentralitySelection]->Fill(phi, eta);
-           fHistRhoVsRCPtRand[fInCentralitySelection]->Fill(pt, RhoVal(phi, fJetRadius, fRho->GetVal())*areaRC);
-           fHistRCPtRand[fInCentralitySelection]->Fill(pt);
-           fHistDeltaPtDeltaPhi2Rand[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt - areaRC*RhoVal(phi, fJetRadius, fRho->GetVal()));
-           fHistDeltaPtDeltaPhi3Rand[fInCentralitySelection]->Fill(PhaseShift(phi-psi3, 3.), pt - areaRC*RhoVal(phi, fJetRadius, fRho->GetVal()));
-       } */
     } 
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillJetHistograms(Double_t vzero[2][2], Double_t* tpc) const
+void AliAnalysisTaskRhoVnModulation::FillJetHistograms(Double_t psi2, Double_t psi3)
 {
     // fill jet histograms
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
@@ -1408,14 +1408,13 @@ void AliAnalysisTaskRhoVnModulation::FillJetHistograms(Double_t vzero[2][2], Dou
         AliEmcalJet* jet = static_cast<AliEmcalJet*>(fJets->At(i));
         if(PassesCuts(jet)) {
             Double_t pt(jet->Pt()), area(jet->Area()), eta(jet->Eta()), phi(jet->Phi());
-            Double_t rho(fLocalRho->GetLocalVal(phi, fJetRadius, fLocalRho->GetVal()));
+            Double_t rho(fLocalRho->GetLocalVal(phi, GetJetContainer()->GetJetRadius(), fLocalRho->GetVal()));
             fHistJetPtRaw[fInCentralitySelection]->Fill(pt);
             fHistJetPt[fInCentralitySelection]->Fill(pt-area*rho);
             if(fFillQAHistograms) fHistJetEtaPhi[fInCentralitySelection]->Fill(eta, phi);
             fHistJetPtArea[fInCentralitySelection]->Fill(pt-area*rho, area);
-            fHistJetPsiTPCPt[fInCentralitySelection]->Fill(PhaseShift(phi-tpc[0], 2.), pt-area*rho);
-            fHistJetPsiVZEROAPt[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[0][0], 2.), pt-area*rho);
-            fHistJetPsiVZEROCPt[fInCentralitySelection]->Fill(PhaseShift(phi-vzero[1][0], 2.), pt-area*rho);
+            fHistJetPsi2Pt[fInCentralitySelection]->Fill(PhaseShift(phi-psi2, 2.), pt-area*rho);
+            fHistJetPsi3Pt[fInCentralitySelection]->Fill(PhaseShift(phi-psi3, 3.), pt-area*rho);
             fHistJetPtConstituents[fInCentralitySelection]->Fill(pt-area*rho, jet->Nch());
             fHistJetEtaRho[fInCentralitySelection]->Fill(eta, pt/area);
             if(fSubtractJetPt) jet->SetPtSub(pt-area*rho);      // if requested, save the subtracted jet pt
@@ -1423,25 +1422,6 @@ void AliAnalysisTaskRhoVnModulation::FillJetHistograms(Double_t vzero[2][2], Dou
     }
 }
 //_____________________________________________________________________________
-void AliAnalysisTaskRhoVnModulation::FillDeltaPhiHistograms(Double_t vzero[2][2], Double_t* tpc) const
-{
-   // fill phi minus psi histograms
-   if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-   if(fTracks) {
-       Int_t iTracks(fTracks->GetEntriesFast());
-       for(Int_t iTPC(0); iTPC < iTracks; iTPC++) {
-           AliVTrack* track = static_cast<AliVTrack*>(fTracks->At(iTPC));
-           if(!PassesCuts(track)) continue;
-           fHistDeltaPhi2VZEROA[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[0][0], 2.));
-           fHistDeltaPhi2VZEROC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[1][0], 2.));
-           fHistDeltaPhi2TPC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-tpc[0], 2.));
-           fHistDeltaPhi3VZEROA[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[0][1], 3.));
-           fHistDeltaPhi3VZEROC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-vzero[1][1], 3.));
-           fHistDeltaPhi3TPC[fInCentralitySelection]->Fill(PhaseShift(track->Phi()-tpc[1], 3.));
-       }
-   }
-}
-//_____________________________________________________________________________
 void AliAnalysisTaskRhoVnModulation::FillQAHistograms(AliVTrack* vtrack) const
 {
     // fill qa histograms for pico tracks
@@ -1481,36 +1461,16 @@ void AliAnalysisTaskRhoVnModulation::FillAnalysisSummaryHistogram() const
 {
     // fill the analysis summary histrogram, saves all relevant analysis settigns
     if(fDebug > 0) printf("__FILE__ = %s \n __LINE __ %i , __FUNC__ %s \n ", __FILE__, __LINE__, __func__);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(1, "fJetRadius"); 
-    fHistAnalysisSummary->SetBinContent(1, fJetRadius);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fPtBiasJetTrack");
-    fHistAnalysisSummary->SetBinContent(2, fPtBiasJetTrack);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fPtBiasJetClus");
-    fHistAnalysisSummary->SetBinContent(3, fPtBiasJetClus);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetPtCut");
-    fHistAnalysisSummary->SetBinContent(4, fJetPtCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetAreaCut");
-    fHistAnalysisSummary->SetBinContent(5, fJetAreaCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fPercAreaCut");
-    fHistAnalysisSummary->SetBinContent(6, fPercAreaCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(7, "fAreaEmcCut");
-    fHistAnalysisSummary->SetBinContent(7, fAreaEmcCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(8, "fJetMinEta");
-    fHistAnalysisSummary->SetBinContent(8, fJetMinEta);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(9, "fJetMaxEta");
-    fHistAnalysisSummary->SetBinContent(9, fJetMaxEta);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(10, "fJetMinPhi");
-    fHistAnalysisSummary->SetBinContent(10, fJetMinPhi);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(11, "fJetMaxPhi");
-    fHistAnalysisSummary->SetBinContent(11, fJetMaxPhi);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(12, "fMaxClusterPt");
-    fHistAnalysisSummary->SetBinContent(12, fMaxClusterPt);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(13, "fMaxTrackPt");
-    fHistAnalysisSummary->SetBinContent(13, fMaxTrackPt);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(14, "fLeadingHadronType");
-    fHistAnalysisSummary->SetBinContent(14, fLeadingHadronType);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(15, "fAnaType");
-    fHistAnalysisSummary->SetBinContent(15, fAnaType);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(2, "fJetRadius");
+    fHistAnalysisSummary->SetBinContent(2, GetJetContainer()->GetJetRadius());
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(3, "fJetEtaMin");
+    fHistAnalysisSummary->SetBinContent(3, GetJetContainer()->GetJetEtaMin());
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(4, "fJetEtaMax");
+    fHistAnalysisSummary->SetBinContent(4, GetJetContainer()->GetJetEtaMax());
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(5, "fJetPhiMin");
+    fHistAnalysisSummary->SetBinContent(5, GetJetContainer()->GetJetPhiMin());
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(6, "fJetPhiMax");
+    fHistAnalysisSummary->SetBinContent(6, GetJetContainer()->GetJetPhiMin());
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(16, "fForceBeamType");
     fHistAnalysisSummary->SetBinContent(16, fForceBeamType);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(17, "fMinCent");
@@ -1523,28 +1483,6 @@ void AliAnalysisTaskRhoVnModulation::FillAnalysisSummaryHistogram() const
     fHistAnalysisSummary->SetBinContent(20, fMaxVz);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(21, "fOffTrigger");
     fHistAnalysisSummary->SetBinContent(21, fOffTrigger);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(22, "fClusPtCut");
-    fHistAnalysisSummary->SetBinContent(22, fClusPtCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(23, "fTrackPtCut");
-    fHistAnalysisSummary->SetBinContent(23, fTrackPtCut);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(24, "fTrackMinEta");
-    fHistAnalysisSummary->SetBinContent(24, fTrackMinEta);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(25, "fTrackMaxEta");
-    fHistAnalysisSummary->SetBinContent(25, fTrackMaxEta);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(26, "fTrackMinPhi");
-    fHistAnalysisSummary->SetBinContent(26, fTrackMinPhi);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(27, "fTrackMaxPhi");
-    fHistAnalysisSummary->SetBinContent(27, fTrackMaxPhi);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(28, "fClusTimeCutLow");
-    fHistAnalysisSummary->SetBinContent(28, fClusTimeCutLow);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(29, "fClusTimeCutUp");
-    fHistAnalysisSummary->SetBinContent(29, fClusTimeCutUp);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(30, "fMinPtTrackInEmcal");
-    fHistAnalysisSummary->SetBinContent(30, fMinPtTrackInEmcal);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(31, "fEventPlaneVsEmcal");
-    fHistAnalysisSummary->SetBinContent(31, fEventPlaneVsEmcal);
-    fHistAnalysisSummary->GetXaxis()->SetBinLabel(32, "fMinEventPlane");
-    fHistAnalysisSummary->SetBinContent(32, fMaxEventPlane);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(33, "fRandomConeRadius");
     fHistAnalysisSummary->SetBinContent(33, fRandomConeRadius);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(34, "fitModulationType");
@@ -1581,6 +1519,10 @@ void AliAnalysisTaskRhoVnModulation::FillAnalysisSummaryHistogram() const
     fHistAnalysisSummary->SetBinContent(49, fSoftTrackMinPt);
     fHistAnalysisSummary->GetXaxis()->SetBinLabel(50, "fSoftTrackMaxPt");
     fHistAnalysisSummary->SetBinContent(50, fSoftTrackMaxPt);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(51, "fMaxCones");
+    fHistAnalysisSummary->SetBinContent(51, fMaxCones);
+    fHistAnalysisSummary->GetXaxis()->SetBinLabel(52, "fUseScaledRho");
+    fHistAnalysisSummary->SetBinContent(52, fUseScaledRho);
 }
 //_____________________________________________________________________________
 void AliAnalysisTaskRhoVnModulation::Terminate(Option_t *)
@@ -1605,6 +1547,13 @@ void AliAnalysisTaskRhoVnModulation::Terminate(Option_t *)
     }
 }
 //_____________________________________________________________________________
+void AliAnalysisTaskRhoVnModulation::SetModulationFit(TF1* fit) 
+{
+// set modulation fit
+    if (fFitModulation) delete fFitModulation;
+    fFitModulation = fit; 
+}
+//_____________________________________________________________________________
 TH1F* AliAnalysisTaskRhoVnModulation::GetResolutionFromOuptutFile(detectorType det, Int_t h, TArrayD* cen)
 {
     // INTERFACE METHOD FOR OUTPUTFILE
@@ -1621,24 +1570,33 @@ TH1F* AliAnalysisTaskRhoVnModulation::GetResolutionFromOuptutFile(detectorType d
         TProfile* temp((TProfile*)fOutputList->FindObject(Form("fProfV%iResolution_%i", h, i)));
         if(!temp) break;
         Double_t a(temp->GetBinContent(3)), b(temp->GetBinContent(5)), c(temp->GetBinContent(7));
+        Double_t d(temp->GetBinContent(9)), e(temp->GetBinContent(10)), f(temp->GetBinContent(11));
         Double_t _a(temp->GetBinError(3)), _b(temp->GetBinError(5)), _c(temp->GetBinError(7));
-        if(a <= 0 || b <= 0 || c <= 0) continue;
+        Double_t _d(temp->GetBinError(9)), _e(temp->GetBinError(10)), _f(temp->GetBinError(11));
+        if(a <= 0 || b <= 0 || c <= 0 || d <= 0 || e <= 0 || f <= 0) continue;
         switch (det) {
             case kVZEROA : {
                 r->SetBinContent(1+i, TMath::Sqrt((a*b)/c));
                 if(i==0) r->SetNameTitle("VZEROA resolution", "VZEROA resolution");
+                r->SetBinError(1+i, TMath::Sqrt(_a*_a+_b*_b+_c*_c));
             } break;
             case kVZEROC : {
                 r->SetBinContent(1+i, TMath::Sqrt((a*c)/b));
                 if(i==0) r->SetNameTitle("VZEROC resolution", "VZEROC resolution");
+                r->SetBinError(1+i, TMath::Sqrt(_a*_a+_b*_b+_c*_c));
             } break;
             case kTPC : {
                 r->SetBinContent(1+i, TMath::Sqrt((b*c)/a));
                 if(i==0) r->SetNameTitle("TPC resolution", "TPC resolution");
+                r->SetBinError(1+i, TMath::Sqrt(_a*_a+_b*_b+_c*_c));
+            } break;
+            case kVZEROComb : {
+                r->SetBinContent(1+i, TMath::Sqrt((d*e)/f));
+                if(i==0) r->SetNameTitle("VZEROComb resolution", "VZEROComb resolution");
+                r->SetBinError(1+i, TMath::Sqrt(_d*_d+_e*_e+_f*_f));
             } break;
             default : break;
         }
-        r->SetBinError(1+i, TMath::Sqrt(_a*_a+_b*_b+_c*_c));
     }
     return r;
 }
@@ -1655,7 +1613,8 @@ TH1F* AliAnalysisTaskRhoVnModulation::CorrectForResolutionDiff(TH1F* v, detector
     Double_t res(1./r->GetBinContent(1+r->FindBin(c)));
     TF1* line = new TF1("line", "pol0", 0, 200);
     line->SetParameter(0, res);
-    return (v->Multiply(line)) ? v : 0x0;
+    v->Multiply(line);
+    return v;
 }
 //_____________________________________________________________________________
 TH1F* AliAnalysisTaskRhoVnModulation::CorrectForResolutionInt(TH1F* v, detectorType det, TArrayD* cen, Int_t h)
@@ -1664,7 +1623,8 @@ TH1F* AliAnalysisTaskRhoVnModulation::CorrectForResolutionInt(TH1F* v, detectorT
     // correct the supplied intetrated vn histogram v for detector resolution
     // integrated vn must have the same centrality binning as the resolotion correction
     TH1F* r(GetResolutionFromOuptutFile(det, h, cen));
-    return (v->Divide(v, r)) ? v : 0x0;
+    v->Divide(v, r);
+    return v;
 }
 //_____________________________________________________________________________
 TH1F* AliAnalysisTaskRhoVnModulation::GetDifferentialQC(TProfile* refCumulants, TProfile* diffCumlants, TArrayD* ptBins, Int_t h)