]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGCF/EBYE/NetParticle/AliAnalysisNetParticleDistribution.cxx
Merge branch 'feature-movesplit'
[u/mrichter/AliRoot.git] / PWGCF / EBYE / NetParticle / AliAnalysisNetParticleDistribution.cxx
index ab36cc5c93347bbd172697540fed88c1d58f38e0..c54431df1bc110341ee3465579f3616eeb4636d3 100644 (file)
@@ -2,31 +2,34 @@
 
 #include "TMath.h"
 #include "TAxis.h"
-#include "TSystem.h" 
 #include "TProfile.h" 
-#include "TH2F.h" 
-#include "TH3F.h" 
-#include "TFile.h" 
-#include "TPRegexp.h"
+#include "TProfile2D.h" 
+#include "TH2D.h" 
+#include "TH3D.h" 
 
 #include "AliStack.h"
 #include "AliMCEvent.h"
 #include "AliMCParticle.h"
+#include "AliESDEvent.h"
 #include "AliESDtrackCuts.h"
-#include "AliESDInputHandler.h"
-#include "AliESDpid.h"
+
 #include "AliCentrality.h"
 #include "AliTracker.h"
-#include "AliAODInputHandler.h"
+
 #include "AliAODEvent.h"
 #include "AliAODTrack.h"
+#include "AliAODMCParticle.h"
 
 #include "AliAnalysisNetParticleDistribution.h"
 
 using namespace std;
 
-// Task for NetParticle checks
-// Author: Jochen Thaeder <jochen@thaeder.de>
+/**
+ * Class for NetParticle Distributions
+ * -- Create input for distributions
+ * Authors: Jochen Thaeder <jochen@thaeder.de>
+ *          Michael Weber <m.weber@cern.ch>
+ */
 
 ClassImp(AliAnalysisNetParticleDistribution)
 
@@ -37,31 +40,19 @@ ClassImp(AliAnalysisNetParticleDistribution)
  */
 
 //________________________________________________________________________
-AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() :
-  fHelper(NULL),
+AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() : 
+  AliAnalysisNetParticleBase("Dist", "Dist"),
 
   fOutList(NULL),
 
-  fESDHandler(NULL),
-  fPIDResponse(NULL),
-  fESD(NULL),
-  fAODHandler(NULL),
-  fAOD(NULL),
-  fIsMC(kFALSE),
-  fMCEvent(NULL),
-  fStack(NULL),
-  fESDTrackCuts(NULL),
-  fEtaMax(0.9),
-  fPtRange(NULL),
-  fAODtrackCutBit(1024),
+  fOrder(8),
+  fNNp(6),
   fNp(NULL),
-  fNCorrNp(2),
-  fCorrNp(NULL),
-  fNMCNp(5),
+  fNpPt(NULL),
+  fNMCNp(7),
   fMCNp(NULL),
-  fNControlMCNp(5),
-  fControlMCNp(NULL),
-
+  fMCNpPt(NULL),
+  fRedFactp(NULL),
   fHnTrackUnCorr(NULL) {
   // Constructor   
   
@@ -72,19 +63,33 @@ AliAnalysisNetParticleDistribution::AliAnalysisNetParticleDistribution() :
 AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() {
   // Destructor
 
-  if (fNp) delete[] fNp;  
+  for (Int_t ii = 0; ii < fNNp; ++ii) 
+    if (fNp[ii]) delete[] fNp[ii];
+  if (fNp) delete[] fNp;
 
-  for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
-    if (fCorrNp[ii]) delete[] fCorrNp[ii];
-  if (fCorrNp) delete[] fCorrNp;
+  for (Int_t ii = 0; ii < 2; ++ii) {
+    for (Int_t kk = 0; kk < 2; ++kk)
+      if (fNpPt[ii][kk]) delete[] fNpPt[ii][kk];      
+    if (fNpPt[ii]) delete[] fNpPt[ii];
+  }
+  if (fNpPt) delete[] fNpPt;
 
+  // -------------------------------------------
   for (Int_t ii = 0; ii < fNMCNp; ++ii) 
     if (fMCNp[ii]) delete[] fMCNp[ii];
   if (fMCNp) delete[] fMCNp;
 
-  for (Int_t ii = 0; ii < fNControlMCNp; ++ii) 
-    if (fControlMCNp[ii]) delete[] fControlMCNp[ii];
-  if (fControlMCNp) delete[] fControlMCNp;
+  for (Int_t ii = 0; ii < 2; ++ii) {
+    for (Int_t kk = 0; kk < 2; ++kk)
+      if (fMCNpPt[ii][kk]) delete[] fMCNpPt[ii][kk];      
+    if (fMCNpPt[ii]) delete[] fMCNpPt[ii];
+  }
+  if (fMCNpPt) delete[] fMCNpPt;
+
+  // -------------------------------------------
+  for (Int_t ii = 0; ii <= fOrder; ++ii) 
+    if (fRedFactp[ii]) delete[] fRedFactp[ii];
+  if (fRedFactp) delete[] fRedFactp;
 
   return;
 }
@@ -96,240 +101,291 @@ AliAnalysisNetParticleDistribution::~AliAnalysisNetParticleDistribution() {
  */
 
 //________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::Initialize(AliAnalysisNetParticleHelper* helper, AliESDtrackCuts* cuts, Bool_t isMC, Float_t *ptRange, Float_t etaMax, Int_t trackCutBit) {
-  // -- Initialize
-  
-  fHelper = helper;
-  fESDTrackCuts = cuts;
-  fIsMC = isMC;
-  fPtRange = ptRange;
-  fEtaMax = etaMax;
-  fAODtrackCutBit = trackCutBit;
-
-  // ------------------------------------------------------------------
-  // -- N particles / N anti-particles
-  // ------------------------------------------------------------------
-  //  Np          : arr[particle]
-  //  CorrNp      : arr[corrSet][particle]
-  //  MCNp        : arr[corrSet][particle] - MC
-  //  ControlMCNp : arr[corrSet][particle] - Control MC
-
-  fNp     = new Float_t[2];
-
-  fCorrNp = new Float_t*[fNCorrNp];
-  for (Int_t ii = 0 ; ii < fNCorrNp; ++ii)
-    fCorrNp[ii] = new Float_t[2];
-
-  fMCNp = new Float_t*[fNMCNp];
-  for (Int_t ii = 0 ; ii < fNMCNp; ++ii)
-    fMCNp[ii] = new Float_t[2];
-
-  fControlMCNp = new Float_t*[fNControlMCNp];
-  for (Int_t ii = 0 ; ii < fNControlMCNp; ++ii)
-    fControlMCNp[ii] = new Float_t[2];
+void AliAnalysisNetParticleDistribution::Process() {
+  // -- Process NetParticle Distributions
 
-  ResetEvent();
+  // -- Fill ESD/AOD tracks
+  ProcessTracks();
+      
+  // -- Fill MC truth particles (missing for AOD MW - However AliVParticle already used)
+  if (fIsMC)
+    ProcessParticles();
 
-  return 0;
+  return;
 }
 
+/*
+ * ---------------------------------------------------------------------------------
+ *                                Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
 //________________________________________________________________________
-void AliAnalysisNetParticleDistribution::CreateHistograms(TList* outList) {
-  // -- Add histograms to outlist
+void AliAnalysisNetParticleDistribution::Init() {
+  // -- Init eventwise
 
-  fOutList = outList;
-  
   // ------------------------------------------------------------------
-  // -- Create net particle histograms
+  // -- N particles / N anti-particles
   // ------------------------------------------------------------------
+  //  Np            : arr[set][particle]
+  //  MCNp          : arr[set][particle] - MC
+  //  Factorials    : arr[order][particle]
 
-  AddHistSet("fHDist", "Uncorrected");
-  AddHistSet("fHDistCorr0", "Corrected [without cross section correction]");
-  AddHistSet("fHDistCorr1", "Corrected [with cross section correction]");
+  //  NpPt          : arr[phiBin][particle][ptBins]
+  //  MCNpPt        : arr[phiBin][particle][ptBins] - MC
+  //  FactorialsPt  : arr[order][particle][ptBins]
   
-  if (fIsMC) {
-    AddHistSet("fHMCrapidity", "MC primary in |y| < 0.5");
-    AddHistSet("fHMCptMin",    "MC primary in |y| + #it{p}_{T} > 0.1");
-    AddHistSet("fHMCpt",       Form("MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
-    AddHistSet("fHMCeta",      Form("MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
-    AddHistSet("fHMCetapt",    Form("MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
-    
-    AddHistSet("fHControlMCLambdarapidity", "Control MC Lambda primary in |y| < 0.5");
-    AddHistSet("fHControlMCLambdaptMin",    "Control MC Lambda primary in |y| + #it{p}_{T} > 0.1");
-    AddHistSet("fHControlMCLambdapt",       Form("Control MC primary in |y| < 0.5 + #it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]));
-    AddHistSet("fHControlMCLambdaeta",      Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f", fEtaMax));
-    AddHistSet("fHControlMCLambdaetapt",    Form("Control MC primary in |y| < 0.5 + |#eta| < %.1f + #it{p}_{T} [%.1f,%.1f]", fEtaMax, fPtRange[0], fPtRange[1]));
+  //   set      ->  fNNp  / fNMCNp
+  //   order    ->  fOrder
+  //   particle ->  2 {0,1}
+  //   phiBin   ->  2 {0,1}
+
+  fNp = new Int_t*[fNNp];
+  for (Int_t ii = 0 ; ii < fNNp; ++ii)
+    fNp[ii] = new Int_t[2];
+
+  fNpPt = new Int_t**[2];
+  for (Int_t ii = 0 ; ii < 2; ++ii) {
+    fNpPt[ii] = new Int_t*[2];
+    for (Int_t kk = 0 ; kk < 2; ++kk)
+      fNpPt[ii][kk] = new Int_t[AliAnalysisNetParticleHelper::fgkfHistNBinsPt];
   }
 
-  // ------------------------------------------------------------------
-  // -- Get Probe Particle Container
-  // ------------------------------------------------------------------
-
-  Double_t dCent[2] = {-0.5, 8.5};
-  Int_t iCent       = 9;
-  
-  Double_t dEta[2]  = {-0.9, 0.9}; // -> 37 bins
-  Int_t iEta        = Int_t((dEta[1]-dEta[0]) / 0.075) +1 ; 
+  // -------------------------------------------
+  fMCNp = new Int_t*[fNMCNp];
+  for (Int_t ii = 0 ; ii < fNMCNp; ++ii)
+    fMCNp[ii] = new Int_t[2];
 
-  Double_t dRap[2]  = {-0.5, 0.5}; 
-  Int_t iRap        = Int_t((dRap[1]-dRap[0]) / 0.075) +1 ; 
+  fMCNpPt = new Int_t**[2];
+  for (Int_t ii = 0 ; ii < 2; ++ii) {
+    fMCNpPt[ii] = new Int_t*[2];
+    for (Int_t kk = 0 ; kk < 2; ++kk)
+      fMCNpPt[ii][kk] = new Int_t[AliAnalysisNetParticleHelper::fgkfHistNBinsPt];
+  }
 
-  Double_t dPhi[2]  = {0.0, TMath::TwoPi()};
-  Int_t iPhi        = 42;
+  // -------------------------------------------
+  fRedFactp = new Double_t*[fOrder+1];
+  for (Int_t ii = 0 ; ii <= fOrder; ++ii)
+    fRedFactp[ii] = new Double_t[2];
+}
 
-  Double_t dPt[2]   = {0.1, 3.0}; 
-  Int_t iPt         = Int_t((dPt[1]-dPt[0]) / 0.075);
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::Reset() {
+  // -- Reset eventwise
 
-  Double_t dSign[2] = {-1.5, 1.5};
-  Int_t iSign       = 3;
+  // -- Reset N particles/anti-particles
+  for (Int_t ii = 0; ii < fNNp; ++ii) 
+    for (Int_t jj = 0; jj < 2; ++jj)
+      fNp[ii][jj] = 0;
 
-  //                          cent:     pt:     sign:     eta:     phi:       y
-  Int_t    binShort[6] = {   iCent,    iPt,    iSign,    iEta,    iPhi,    iRap};
-  Double_t minShort[6] = {dCent[0], dPt[0], dSign[0], dEta[0], dPhi[0], dRap[0]};
-  Double_t maxShort[6] = {dCent[1], dPt[1], dSign[1], dEta[1], dPhi[1], dRap[1]};
+  for (Int_t ii = 0; ii < 2; ++ii) 
+    for (Int_t jj = 0; jj < 2; ++jj)
+      for (Int_t kk = 0; kk < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++kk)
+       fNpPt[ii][jj][kk] = 0;
 
-  // -- UnCorrected
-  fOutList->Add(new THnSparseF("fHnTrackUnCorr", "cent:pt:sign:eta:phi:y", 6, binShort, minShort, maxShort));  
-  fHnTrackUnCorr = static_cast<THnSparseF*>(fOutList->Last());
-  fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality");
-  fHnTrackUnCorr->GetAxis(1)->SetTitle("#it{p}_{T} (GeV/#it{c})");
-  fHnTrackUnCorr->GetAxis(2)->SetTitle("sign");
-  fHnTrackUnCorr->GetAxis(3)->SetTitle("#eta");
-  fHnTrackUnCorr->GetAxis(4)->SetTitle("#varphi");
-  fHnTrackUnCorr->GetAxis(5)->SetTitle("#it{y}");
+  // -- Reset N MC particles/anti-particles
+  for (Int_t ii = 0; ii < fNMCNp; ++ii) 
+    for (Int_t jj = 0; jj < 2; ++jj)
+      fMCNp[ii][jj] = 0;
 
-  fHelper->BinLogAxis(fHnTrackUnCorr, 1);
+  for (Int_t ii = 0; ii < 2; ++ii) 
+    for (Int_t jj = 0; jj < 2; ++jj)
+      for (Int_t kk = 0; kk < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++kk)
+       fMCNpPt[ii][jj][kk] = 0;
 
-  // ------------------------------------------------------------------
+  // -- Reset reduced factorials for particles/anti-particles
+  for (Int_t ii = 0; ii <= fOrder; ++ii) 
+    for (Int_t jj = 0; jj < 2; ++jj)
+      fRedFactp[ii][jj] = 1.;
 
-  return;
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::SetupEvent(AliESDInputHandler *esdHandler, AliAODInputHandler *aodHandler,  AliMCEvent *mcEvent) {
-  // -- Setup Event
+void AliAnalysisNetParticleDistribution::CreateHistograms() {
+  // -- Add histograms to outlist
 
-  ResetEvent();
+  // ------------------------------------------------------------------
+  // -- Get Probe Particle Container
+  // ------------------------------------------------------------------
 
-  // -- Get ESD objects
-  if(esdHandler){
-    fESDHandler  = esdHandler;
-    fPIDResponse = esdHandler->GetPIDResponse();
-    fESD         = fESDHandler->GetEvent();
-  }
+  Int_t    binHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistNBinsCent, AliAnalysisNetParticleHelper::fgkfHistNBinsEta,          // cent |  eta
+                            AliAnalysisNetParticleHelper::fgkfHistNBinsRap,  AliAnalysisNetParticleHelper::fgkfHistNBinsPhi,          //    y |  phi
+                            AliAnalysisNetParticleHelper::fgkfHistNBinsPt,   AliAnalysisNetParticleHelper::fgkfHistNBinsSign};        //   pt | sign
+  
+  Double_t minHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeEta[0], 
+                            AliAnalysisNetParticleHelper::fgkfHistRangeRap[0],  AliAnalysisNetParticleHelper::fgkfHistRangePhi[0], 
+                            AliAnalysisNetParticleHelper::fgkfHistRangePt[0],   AliAnalysisNetParticleHelper::fgkfHistRangeSign[0]};
+  
+  Double_t maxHnUnCorr[6] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[1], AliAnalysisNetParticleHelper::fgkfHistRangeEta[1], 
+                            AliAnalysisNetParticleHelper::fgkfHistRangeRap[1],  AliAnalysisNetParticleHelper::fgkfHistRangePhi[1], 
+                            AliAnalysisNetParticleHelper::fgkfHistRangePt[1],   AliAnalysisNetParticleHelper::fgkfHistRangeSign[1]};
+  
+  // -- UnCorrected
+  fOutList->Add(new THnSparseD("hnTrackUnCorr", "cent:eta:y:phi:pt:sign", 6, binHnUnCorr, minHnUnCorr, maxHnUnCorr));  
+  fHnTrackUnCorr = static_cast<THnSparseD*>(fOutList->Last());
+  fHnTrackUnCorr->Sumw2(); 
+  fHnTrackUnCorr->GetAxis(0)->SetTitle("centrality");
+  fHnTrackUnCorr->GetAxis(1)->SetTitle("#eta");
+  fHnTrackUnCorr->GetAxis(2)->SetTitle("#it{y}");
+  fHnTrackUnCorr->GetAxis(3)->SetTitle("#varphi");
+  fHnTrackUnCorr->GetAxis(4)->SetTitle("#it{p}_{T} (GeV/#it{c})");
+  fHnTrackUnCorr->GetAxis(5)->SetTitle("sign");
 
-  // -- Get AOD objects
-  else if(aodHandler){
-    fAODHandler  = aodHandler;
-    fPIDResponse = aodHandler->GetPIDResponse();
-    fAOD         = fAODHandler->GetEvent();
-  }
+  fHelper->BinLogAxis(fHnTrackUnCorr, 4, fESDTrackCuts);
 
-  // -- Get MC objects
-  fMCEvent     = mcEvent;
-  if (fMCEvent)
-    fStack     = fMCEvent->Stack();
+  for (Int_t idx = 1 ; idx <= fHnTrackUnCorr->GetAxis(4)->GetNbins(); ++idx)
+    printf("%02d |  %f > %f < %f\n", idx, fHnTrackUnCorr->GetAxis(4)->GetBinLowEdge(idx), fHnTrackUnCorr->GetAxis(4)->GetBinCenter(idx), fHnTrackUnCorr->GetAxis(4)->GetBinUpEdge(idx));
 
-  return 0;
-}
+  // ------------------------------------------------------------------
+  // -- Create net particle histograms
+  // ------------------------------------------------------------------
 
-//________________________________________________________________________
-void AliAnalysisNetParticleDistribution::ResetEvent() {
-  // -- Reset event
-  
-  // -- Reset ESD Event
-  fESD       = NULL;
+  // -- Get ranges for pt and eta 
+  Float_t etaRange[2];
+  fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 
-  // -- Reset AOD Event
-  fAOD       = NULL;
+  Float_t ptRange[2];
+  fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
 
-  // -- Reset MC Event
-  if (fIsMC)
-    fMCEvent = NULL;
+  // ------------------------------------------------------------------
 
-  // -- Reset N particles/anti-particles
-  for (Int_t jj = 0; jj < 2; ++jj)
-    fNp[jj] = 0.;
+  TString sTitle("");
+  sTitle = (fHelper->GetUsePID()) ? Form("|y|<%.1f", fHelper->GetRapidityMax()) : Form("|#eta|<%.1f", etaRange[1]);
+
+  // -- centrality dependent
+  AddHistSetCent("Dist",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));
+  AddHistSetCent("DistTPC",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
+  AddHistSetCent("DistTOF",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
+
+#if USE_PHI
+  AddHistSetCent("Distphi",    Form("%s,#it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                    sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+  AddHistSetCent("DistTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                    sTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+  AddHistSetCent("DistTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                    sTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
+
+  // -- centrality + pt dependent
+  AddHistSetCentPt("PtDist",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sTitle.Data(), ptRange[0], ptRange[1]));  
+#if USE_PHI
+  AddHistSetCentPt("PtDistphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                     sTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
   
-  // -- Reset N corrected particles/anti-particles
-  for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
-    for (Int_t jj = 0; jj < 2; ++jj)
-      fCorrNp[ii][jj] = 0.;
+  if (fIsMC) {
+    TString sMCTitle("");
+    sMCTitle = (fHelper->GetUsePID()) ?  Form("MC primary in |y| < %.1f", fHelper->GetRapidityMax()) : Form("MC primary in |#eta| < %.1f", etaRange[1]);
+
+    // -- centrality dependent
+    AddHistSetCent("MC",         Form("%s", sTitle.Data()));
+
+    AddHistSetCent("MCpt",       Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1]));
+    AddHistSetCent("MCptTPC",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired()));
+    AddHistSetCent("MCptTOF",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1]));
+
+#if USE_PHI
+    AddHistSetCent("MCptphi",    Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                      sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+    AddHistSetCent("MCptTPCphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                      sMCTitle.Data(), ptRange[0], fHelper->GetMinPtForTOFRequired(), fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+    AddHistSetCent("MCptTOFphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                      sMCTitle.Data(), fHelper->GetMinPtForTOFRequired(), ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
+
+    // -- centrality + pt dependent
+    AddHistSetCentPt("PtMCpt",    Form("%s, #it{p}_{T} [%.1f,%.1f]", sMCTitle.Data(), ptRange[0], ptRange[1]));
+#if USE_PHI
+    AddHistSetCentPt("PtMCptphi", Form("%s, #it{p}_{T} [%.1f,%.1f], #varphi [%.2f,%.2f]", 
+                                       sMCTitle.Data(), ptRange[0], ptRange[1], fHelper->GetPhiMin(), fHelper->GetPhiMax()));
+#endif
+  }
 
-  // -- Reset N MC particles/anti-particles
-  for (Int_t ii = 0; ii < fNMCNp; ++ii) 
-    for (Int_t jj = 0; jj < 2; ++jj)
-      fMCNp[ii][jj] = 0.;
+  // ------------------------------------------------------------------
 
-  // -- Reset N control MC particles/anti-particles
-  for (Int_t ii = 0; ii < fNControlMCNp; ++ii) 
-    for (Int_t jj = 0; jj < 2; ++jj)
-      fControlMCNp[ii][jj] = 0.;
+  return;
 }
 
+/*
+ * ---------------------------------------------------------------------------------
+ *                                Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
 //________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::Process() {
-  // -- Process NetParticle Distributions
+Int_t AliAnalysisNetParticleDistribution::ProcessTracks() {
+  // -- Process ESD/AOD tracks and fill histograms
 
-  // -- Fill ESD tracks
-  if (fESD) ProcessESDTracks();
+  // -- Get ranges for AOD particles
+  Float_t etaRange[2];
+  fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
 
-  // -- Fill AOD tracks
-  else if (fAOD) ProcessAODTracks();
+  Float_t ptRange[2];
+  fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+  
+  // -- Track Loop
+  for (Int_t idxTrack = 0; idxTrack < fNTracks; ++idxTrack) {
     
-  // -- Fill MC truth particles (missing for AOD XXX)
-  if (fIsMC)  {
-    ProcessStackParticles();
-    ProcessStackControlParticles();
-  }
-
-  return 0;
-}
-
-//________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() {
-  // -- Process ESD tracks and fill histograms
-
-  for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
-    AliESDtrack *track = fESD->GetTrack(idxTrack); 
+    AliVTrack *track = (fESD) ? static_cast<AliVTrack*>(fESD->GetTrack(idxTrack)) : static_cast<AliVTrack*>(fAOD->GetTrack(idxTrack)); 
 
     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
     // -- Check track cuts
     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-    
-    // -- Check if charged track is accepted for basic parameters
+
+    // -- Check if track is accepted for basic parameters
     if (!fHelper->IsTrackAcceptedBasicCharged(track))
       continue;
     
-    // -- Check if accepted
-    if (!fESDTrackCuts->AcceptTrack(track)) 
+    // -- Check if accepted - ESD
+    if (fESD && !fESDTrackCuts->AcceptTrack(dynamic_cast<AliESDtrack*>(track)))
       continue;
-
-    // -- Check if accepted in rapidity window
+    
+    // -- Check if accepted - AOD
+    if (fAOD){
+      AliAODTrack * trackAOD = dynamic_cast<AliAODTrack*>(track);
+      
+      if (!trackAOD) {
+       AliError("Pointer to dynamic_cast<AliAODTrack*>(track) = ZERO");
+       continue;
+      }
+      if (!trackAOD->TestFilterBit(fAODtrackCutBit))
+       continue;
+      
+      // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
+      if(!(track->Pt() > ptRange[0] && track->Pt() <= ptRange[1] && TMath::Abs(track->Eta()) <= etaRange[1]))
+       continue;
+    }
+    
+    // -- Check if accepted in rapidity window -- for identified particles
+    //    for charged - eta check is done in the trackcuts
     Double_t yP;
-    if (!fHelper->IsTrackAcceptedRapidity(track, yP))
+    if (fHelper->GetUsePID() && !fHelper->IsTrackAcceptedRapidity(track, yP))
       continue;
 
-    // -- Check if accepted bt PID from TPC or TPC+TOF
-    Double_t pid[2];
-    if (!fHelper->IsTrackAcceptedPID(track, pid))
+    // -- Check if accepted with thighter DCA cuts
+    // -- returns kTRUE for AODs for now : MW
+    if (!fHelper->IsTrackAcceptedDCA(track))
       continue;
 
-    // -- Check if accepted with thighter DCA cuts
-    if (fHelper->IsTrackAcceptedDCA(track))
+    // -- Check if accepted by PID from TPC or TPC+TOF
+    Double_t pid[3];
+    if (!fHelper->IsTrackAcceptedPID(track, pid))
       continue;
-    
+
     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
     // -- Fill Probe Particle
     // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 
+    if (!fHelper->GetUsePID())  
+      yP = track->Eta();
+    
     Double_t aTrack[6] = {
-      Double_t(fHelper->GetCentralityBin()),       //  0 centrality 
-      track->Pt(),                    //  1 pt
-      track->GetSign(),               //  2 sign
-      track->Eta(),                   //  3 eta
-      track->Phi(),                   //  4 phi
-      yP                              //  5 rapidity
+      Double_t(fCentralityBin),               //  0 centrality 
+      track->Eta(),                           //  1 eta
+      yP,                                     //  2 rapidity
+      track->Phi(),                           //  3 phi
+      track->Pt(),                            //  4 pt
+      static_cast<Double_t>(track->Charge())                         //  5 sign
     };
     
     fHnTrackUnCorr->Fill(aTrack);
@@ -338,116 +394,85 @@ Int_t AliAnalysisNetParticleDistribution::ProcessESDTracks() {
     // ------------------------------------------------------------------
     //  idxPart = 0 -> anti particle
     //  idxPart = 1 -> particle
+    Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
 
-    Int_t idxPart = (track->GetSign() < 0) ? 0 : 1;
-    fNp[idxPart] += 1.;
-    
-    for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
-      fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);      
-    
-  } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
-
-  // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-  // -- Fill Particle Fluctuation Histograms
-  // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+    //  idxPhi  = 0 -> Full Acceptance
+    //  idxPhi  = 1 -> Partial Acceptance
+    Int_t idxPhi = 0;
 
-  // -- Uncorrected
-  FillHistSet("fHDist", fNp);
-    
-  // -- Corrected 
-  for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) 
-    FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]); 
+    //  idxPt       -> [0, nBins-1]
+    Int_t idxPt = static_cast<TAxis*>(fHnTrackUnCorr->GetAxis(4))->FindBin(track->Pt())-1;
 
-  return 0;
-}
+    // -- using pt Bins
+    fNpPt[idxPhi][idxPart][idxPt] += 1;
 
-//________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::ProcessAODTracks() {
-  // -- Process ESD tracks and fill histograms
+    // -- in pt Range
+    fNp[0][idxPart] += 1;
 
-  for (Int_t idxTrack = 0; idxTrack < fAOD->GetNumberOfTracks(); ++idxTrack) {
-    AliAODTrack *track = (AliAODTrack*)fAOD->GetTrack(idxTrack); 
+    // -- in TPC pt Range
+    if (track->Pt() <= fHelper->GetMinPtForTOFRequired())
+      fNp[1][idxPart] += 1;
 
-    // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-    // -- Check track cuts
-    // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+    // -- in TPC+TOF pt Range
+    if (track->Pt() > fHelper->GetMinPtForTOFRequired())
+      fNp[2][idxPart] += 1;
 
-    // -- Check if charged track is accepted for basic parameters
-    if (!fHelper->IsTrackAcceptedBasicCharged(track))
+    // -- check phi range ----------------------------------------------------------
+#if USE_PHI
+    if(!fHelper->IsTrackAcceptedPhi(track))
       continue;
+    idxPhi = 1;
 
-    // -- Check if accepted
-    if(!track->TestFilterBit(fAODtrackCutBit)) 
-      continue;
-    
-    // -- Check if in pT and eta range (is done in ESDTrackCuts for ESDs)
-    if(!(track->Pt() > fPtRange[0] && track->Pt() < fPtRange[1] && TMath::Abs(track->Eta()) <= fEtaMax))
-      continue;
+    // -- using pt Bins
+    fNpPt[idxPhi][idxPart][idxPt] += 1;
 
-    // -- Check if accepted in rapidity window
-    Double_t yP;
-    if (!fHelper->IsTrackAcceptedRapidity(track, yP))
-      continue;
-
-    // -- Check if accepted bt PID from TPC or TPC+TOF
-    Double_t pid[2];
-    if (!fHelper->IsTrackAcceptedPID(track, pid))
-      continue;
-
-    // -- Check if accepted with thighter DCA cuts XXX
-    // if (fHelper->IsTrackAcceptedDCA(track))
-    //  continue;
-    
-    // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-    // -- Fill Probe Particle
-    // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-
-    Double_t aTrack[6] = {
-      Double_t(fHelper->GetCentralityBin()),       //  0 centrality 
-      track->Pt(),                    //  1 pt
-      track->Charge(),                //  2 sign
-      track->Eta(),                   //  3 eta
-      track->Phi(),                   //  4 phi
-      yP                              //  5 rapidity
-    };
-    
-    fHnTrackUnCorr->Fill(aTrack);
-    
-    // -- Count particle / anti-particle 
-    // ------------------------------------------------------------------
-    //  idxPart = 0 -> anti particle
-    //  idxPart = 1 -> particle
+    // -- in pt Range
+    fNp[3][idxPart] += 1;
 
-    Int_t idxPart = (track->Charge() < 0) ? 0 : 1;
-    fNp[idxPart] += 1.;
-    
-    for (Int_t ii = 0; ii < fNCorrNp; ++ii) 
-      fCorrNp[ii][idxPart] += fHelper->GetTrackbyTrackCorrectionFactor(aTrack, ii);      
+    // -- in TPC pt Range
+    if (track->Pt() <= fHelper->GetMinPtForTOFRequired()) 
+      fNp[4][idxPart] += 1;
     
+    // -- in TPC+TOF pt Range
+    if (track->Pt() > fHelper->GetMinPtForTOFRequired())
+      fNp[5][idxPart] += 1;
+#endif
   } // for (Int_t idxTrack = 0; idxTrack < fESD->GetNumberOfTracks(); ++idxTrack) {
 
   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
   // -- Fill Particle Fluctuation Histograms
   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 
-  // -- Uncorrected
-  FillHistSet("fHDist", fNp);
-    
-  // -- Corrected 
-  for (Int_t ii = 0 ; ii < fNCorrNp; ++ii) 
-    FillHistSet(Form("fHDistCorr%d", ii), fCorrNp[ii]); 
+  FillHistSetCent("Dist",        0, kFALSE);
+  FillHistSetCent("DistTPC",     1, kFALSE);
+  FillHistSetCent("DistTOF",     2, kFALSE);
+
+  FillHistSetCentPt("PtDist",    0, kFALSE);
 
+#if USE_PHI
+  FillHistSetCent("Distphi",     3, kFALSE);
+  FillHistSetCent("DistTPCphi",  4, kFALSE);
+  FillHistSetCent("DistTOFphi",  5, kFALSE);
+
+  FillHistSetCentPt("PtDistphi", 1, kFALSE);
+#endif
   return 0;
 }
 
 //________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() {
+Int_t AliAnalysisNetParticleDistribution::ProcessParticles() {
   // -- Process primary particles from the stack and fill histograms
-
-  Int_t pdgCode    = AliPID::ParticleCode(fHelper->GetParticleSpecies());
   
+  Float_t etaRange[2];
+  fESDTrackCuts->GetEtaRange(etaRange[0],etaRange[1]);
+
+  Float_t ptRange[2];
+  fESDTrackCuts->GetPtRange(ptRange[0],ptRange[1]);
+
   for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
-    TParticle* particle = fStack->Particle(idxMC);
+    AliVParticle* particle = (fESD) ? fMCEvent->GetTrack(idxMC) : NULL;
+
     if (!particle) 
       continue;
 
@@ -456,229 +481,514 @@ Int_t AliAnalysisNetParticleDistribution::ProcessStackParticles() {
       continue;
     
     // -- Check if particle / anti-particle
-    if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
+    if (fHelper->GetUsePID() && TMath::Abs(particle->PdgCode()) != fPdgCode)
       continue;
     
-    // -- Get particle : 0 anti-particle / 1 particle
-    Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
-
-    // >> NOW only anti-particle / particle 
-    // >> With idxPart
-    
-    // -- Check rapidity window
+    // -- Check rapidity window -- for identfied particles
     Double_t yMC;
-    if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
+    if (fHelper->GetUsePID() && !fHelper->IsParticleAcceptedRapidity(particle, yMC))
+      continue;
+
+    // -- Check eta window -- for charged particles
+    if (!fHelper->GetUsePID() && TMath::Abs(particle->Eta()) > etaRange[1])
       continue;
     
-    fMCNp[0][idxPart] += 1.;         // -> MCrapidity
-    
-    // -- Check acceptance
-    if (particle->Pt() > 0.1 )
-      fMCNp[1][idxPart] += 1.;       // -> MCptMin
-    
-    if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
-      fMCNp[2][idxPart] += 1.;       // -> MCpt
-    
-    if (TMath::Abs(particle->Eta()) <= fEtaMax)
-      fMCNp[3][idxPart] += 1.;       // -> MCeta
+    // -- Count particle / anti-particle 
+    // ------------------------------------------------------------------
+    //  idxPart = 0 -> anti particle
+    //  idxPart = 1 -> particle
+    Int_t idxPart = (particle->PdgCode() < 0) ? 0 : 1;
+
+    //  idxPhi  = 0 -> Full Acceptance
+    //  idxPhi  = 1 -> Partial Acceptance
+    Int_t idxPhi = 0;
+
+    //  idxPt       -> [0, nBins-1]
+    Int_t idxPt = static_cast<TAxis*>(fHnTrackUnCorr->GetAxis(4))->FindBin(particle->Pt()) - 1;
+
+    // -- MCrapidity for identfied particles
+    //    MCeta for charged particles
+    fMCNp[0][idxPart] += 1.;        
+     
+    // -- Check main pt window
+    if (!(particle->Pt() > ptRange[0] && particle->Pt() <= ptRange[1]))
+      continue;
+
+    // -- using pt Bin
+    fMCNpPt[idxPhi][idxPart][idxPt] += 1;
+
+    // -- in pt Range
+    fMCNp[1][idxPart] += 1;
     
-    if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
-      fMCNp[4][idxPart] += 1.;       // -> MCetapt
+    // -- in TPC pt Range
+    if (particle->Pt() <= fHelper->GetMinPtForTOFRequired())
+      fMCNp[2][idxPart] += 1;
+
+    // -- in TPC+TOF pt Range
+    if (particle->Pt() > fHelper->GetMinPtForTOFRequired())
+      fMCNp[3][idxPart] += 1;
+
+    // -- check phi range ----------------------------------------------------------
+#if USE_PHI
+    if(!fHelper->IsParticleAcceptedPhi(particle))
+      continue;
+    idxPhi = 1;
+
+    // -- using pt Bin
+    fMCNpPt[idxPhi][idxPart][idxPt] += 1;
     
+    // -- in pt Range
+    fMCNp[4][idxPart] += 1;
+
+    // -- in TPC pt Range
+    if (particle->Pt() <= fHelper->GetMinPtForTOFRequired())
+      fMCNp[5][idxPart] += 1;
+
+    // -- in TPC+TOF pt Range
+    if (particle->Pt() > fHelper->GetMinPtForTOFRequired())
+      fMCNp[6][idxPart] += 1;
+#endif
   } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
   
   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
   // -- Fill Particle Fluctuation Histograms
   // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
 
-  FillHistSet("fHMCrapidity", fMCNp[0]);
-  FillHistSet("fHMCptMin",    fMCNp[1]);
-  FillHistSet("fHMCpt",       fMCNp[2]);
-  FillHistSet("fHMCeta",      fMCNp[3]);
-  FillHistSet("fHMCetapt",    fMCNp[4]);
+  FillHistSetCent("MC",          0, kTRUE);
+  FillHistSetCent("MCpt",        1, kTRUE);
+  FillHistSetCent("MCptTPC",     2, kTRUE);
+  FillHistSetCent("MCptTOF",     3, kTRUE);
+
+  FillHistSetCentPt("PtMCpt",    0, kTRUE);
+
+#if USE_PHI
+  FillHistSetCent("MCptphi",     4, kTRUE);
+  FillHistSetCent("MCptTPCphi",  5, kTRUE);
+  FillHistSetCent("MCptTOFphi",  6, kTRUE);
+
+  FillHistSetCentPt("PtMCptphi", 1, kTRUE);
+#endif
 
   return 0;
 }
 
+/*
+ * ---------------------------------------------------------------------------------
+ *                            Helper Methods - private
+ * ---------------------------------------------------------------------------------
+ */
+
 //________________________________________________________________________
-Int_t AliAnalysisNetParticleDistribution::ProcessStackControlParticles() {
-  // -- Process primary particles from the stack and fill histograms
+void  AliAnalysisNetParticleDistribution::AddHistSetCent(const Char_t *name, const Char_t *title)  {
+  // -- Add histogram sets for particle and anti-particle
+  //    dependence : centrality
 
-  Int_t pdgCode    = fHelper->GetControlParticleCode();
-  Bool_t isNeutral  = fHelper->IsControlParticleNeutral();
-  const Char_t* name = fHelper->GetControlParticleName().Data();
+  TString sName(name);
+  TString sTitle(title);
 
-  for (Int_t idxMC = 0; idxMC < fStack->GetNprimary(); ++idxMC) {
-    TParticle* particle = fStack->Particle(idxMC);
-    if (!particle) 
-      continue;
-    
-    // -- Check basic MC properties -> neutral or charged physical primary
-    if (isNeutral) {
-      if (!fHelper->IsParticleAcceptedBasicNeutral(particle, idxMC))
-       continue;
-    }
-    else {
-      if (!fHelper->IsParticleAcceptedBasicCharged(particle, idxMC))
-       continue;
-    }
-    
-    // -- Check if particle / anti-particle
-    if (TMath::Abs(particle->GetPdgCode()) != pdgCode)
-      continue;
+  // -- Add List
+  fOutList->Add(new TList);
+  TList *list = static_cast<TList*>(fOutList->Last());
+  list->SetName(Form("f%s", name));
+  list->SetOwner(kTRUE);
 
-    // -- Get particle : 0 anti-particle / 1 particle
-    Int_t idxPart = (particle->GetPdgCode() < 0) ? 0 : 1;
+  // -- Get Bin Ranges
+  Int_t nBinsCent         =  AliAnalysisNetParticleHelper::fgkfHistNBinsCent;
+  Double_t centBinRange[] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeCent[1]};
 
-    // >> NOW only anti-particle / particle 
-    // >> With idxPart
-    
-    // -- Check rapidity window
-    Double_t yMC;
-    if (!fHelper->IsParticleAcceptedRapidity(particle, yMC))
-      continue;
+  // -- Create Titles
+  TString sNetTitle(Form("N_{%s} - N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data()));
+  TString sSumTitle(Form("N_{%s} + N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data()));
 
-    fControlMCNp[0][idxPart] += 1.;         // -> ControlMCrapidity
+  // -- Add Particle / Anti-Particle Distributions
+  for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+    list->Add(new TH2D(Form("h%s%s", name, fHelper->GetParticleName(idxPart).Data()), 
+                      Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()),
+                      nBinsCent, centBinRange[0], centBinRange[1], 2601, -0.5, 2600.49));
+  } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+   
+  // -- Add NetParticle Distributions
+  list->Add(new TH2D(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()), 
+                    Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()), 
+                      nBinsCent, centBinRange[0], centBinRange[1], 601, -300.5, 300.49));
 
-    // -- Check acceptance
-    if (particle->Pt() > 0.1 )
-      fControlMCNp[1][idxPart] += 1.;       // -> ControlMCptMin
-    
-    if (particle->Pt() > fPtRange[0] && particle->Pt() <= fPtRange[1])
-      fControlMCNp[2][idxPart] += 1.;       // -> ControlMCpt
-    
-    if (TMath::Abs(particle->Eta()) <= fEtaMax)
-      fControlMCNp[3][idxPart] += 1.;       // -> ControlMCeta
-    
-    if (particle->Pt() > fPtRange[0] && particle->Pt() > fPtRange[1] && TMath::Abs(particle->Eta()) <= fEtaMax)
-      fControlMCNp[4][idxPart] += 1.;       // -> ControlMCetapt
-  } // for (Int_t idxMC = 0; idxMC < nPart; ++idxMC) {
+  // -- Add NetParticle vs SumParticle
+  list->Add(new TH2D(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), 
+                    Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), 
+                      nBinsCent, centBinRange[0], centBinRange[1], 41, -2.5, 2.49));
 
-  // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
-  // -- Fill Particle Fluctuation Histograms
-  // -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- -- --
+  // -----------------------------------------------------------------------------------------------
 
-  FillHistSet(Form("fHControlMC%srapidity", name), fControlMCNp[0], 0);
-  FillHistSet(Form("fHControlMC%sptMin",    name), fControlMCNp[1], 1);
-  FillHistSet(Form("fHControlMC%spt",       name), fControlMCNp[2], 2);
-  FillHistSet(Form("fHControlMC%seta",      name), fControlMCNp[3], 3);
-  FillHistSet(Form("fHControlMC%setapt",    name), fControlMCNp[4], 4);
+  // -- Add Particle / Anti-Particle Distributions
+  for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+    list->Add(new TH2D(Form("h%s%sX", name, fHelper->GetParticleName(idxPart).Data()), 
+                      Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()),
+                      nBinsCent, centBinRange[0], centBinRange[1], 2601, -0.5, 2600.49));
+  } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+   
+  // -- Add NetParticle Distributions
+  list->Add(new TH2D(Form("h%sNet%sX", name, fHelper->GetParticleName(1).Data()), 
+                    Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()), 
+                      nBinsCent, centBinRange[0], centBinRange[1], 601, -300.5, 300.49));
+
+  // -- Add NetParticle vs SumParticle
+  list->Add(new TH2D(Form("h%sNet%sOverSumX", name, fHelper->GetParticleName(1).Data()), 
+                    Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), 
+                      nBinsCent, centBinRange[0], centBinRange[1], 41, -2.5, 2.49));
+
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <NetParticle^k>
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idx = 1; idx <= fOrder; ++idx) {
+    list->Add(new TProfile(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idx), 
+                          Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+                          nBinsCent, centBinRange[0], centBinRange[1]));
+  }
 
-  return 0;
-}
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <NetParticle^k> for every SubSample
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) {
+    for (Int_t idx = 1; idx <= fOrder; ++idx) {
+      list->Add(new TProfile(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idx, idxSub), 
+                            Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+                            nBinsCent, centBinRange[0], centBinRange[1]));
+    }
+  }
+  
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <f_ik>
+  // -----------------------------------------------------------------------------------------------
+  list->Add(new TList);
+  TList *fikList = static_cast<TList*>(list->Last());
+  fikList->SetName(Form("f%sFik",name));
+  fikList->SetOwner(kTRUE);
+
+  for (Int_t ii = 0; ii <= fOrder; ++ii) {
+    for (Int_t kk = 0; kk <= fOrder; ++kk) {
+      fikList->Add(new TProfile(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk),
+                               Form("f_{%02d%02d} : %s;Centrality;f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
+                               nBinsCent, centBinRange[0], centBinRange[1]));
+    }
+  }
 
-////////////////////////////////////////////////////////////////////////////////////////////
+  // -- Add counter for number of Non-zero entries
+  for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) 
+    fikList->Add(new TH2D(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), idxCent),
+                         Form("f_ik counts : %s Cent %d;i;k;", sTitle.Data(), idxCent),
+                         fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49));
+  
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <f_ik> for every SubSample
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) {
+    list->Add(new TList);
+    TList *fikListSub = static_cast<TList*>(list->Last());
+    fikListSub->SetName(Form("f%sFik_%02d",name, idxSub));
+    fikListSub->SetOwner(kTRUE);
+
+    for (Int_t ii = 0; ii <= fOrder; ++ii) {
+      for (Int_t kk = 0; kk <= fOrder; ++kk) {
+       fikListSub->Add(new TProfile(Form("p%sNet%sF%02d%02d_%02d", name, fHelper->GetParticleName(1).Data(), ii, kk, idxSub),
+                                    Form("f_{%02d%02d} : %s;Centrality;f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk),
+                                    nBinsCent, centBinRange[0], centBinRange[1]));
+      }
+    }
 
+    // -- Add counter for number of Non-zero entries
+    for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) 
+      fikListSub->Add(new TH2D(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), idxCent, idxSub),
+                              Form("f_ik counts : %s Cent %d;i;k;", sTitle.Data(), idxCent),
+                              fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49));
+  }
+  
+  return;
+}
 
 //________________________________________________________________________
-void  AliAnalysisNetParticleDistribution::AddHistSet(const Char_t *name, const Char_t *title)  {
+void  AliAnalysisNetParticleDistribution::AddHistSetCentPt(const Char_t *name, const Char_t *title)  {
   // -- Add histogram sets for particle and anti-particle
+  //    dependence : centrality and pt 
+  
+  TString sName(name);
+  TString sTitle(title);
 
-  Int_t pid = fHelper->GetParticleSpecies();
-  if( pid < 0 || pid > AliPID::kSPECIES){
-    AliError("Particle ID not in AliPID::kSPECIES");
-    return;
-  }
-
+  // -- Add List
   fOutList->Add(new TList);
   TList *list = static_cast<TList*>(fOutList->Last());
-  list->SetName(name) ;
+  list->SetName(Form("f%s", name));
   list->SetOwner(kTRUE);
 
-  TString sName(name);
+  // -- Get Bin Ranges
+  Int_t nBinsPt           =  AliAnalysisNetParticleHelper::fgkfHistNBinsPt;
+  Double_t ptBinRange[]   = {-0.5, AliAnalysisNetParticleHelper::fgkfHistNBinsPt - 0.5};
 
-  TString sPtTitle("");
-  if (!sName.Contains("fHMC") || !sName.Contains("fHControlMC"))
-    sPtTitle += Form("#it{p}_{T} [%.1f,%.1f]", fPtRange[0], fPtRange[1]);
-  
+  Int_t nBinsCent         =  AliAnalysisNetParticleHelper::fgkfHistNBinsCent;
+  Double_t centBinRange[] = {AliAnalysisNetParticleHelper::fgkfHistRangeCent[0], AliAnalysisNetParticleHelper::fgkfHistRangeCent[1]};
+
+  // -- Create Titles
+  TString sNetTitle(Form("N_{%s} - N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data()));
+  TString sSumTitle(Form("N_{%s} + N_{%s}", fHelper->GetParticleTitleLatex(1).Data(), fHelper->GetParticleTitleLatex(0).Data()));
+
+  // -- Add Particle / Anti-Particle Distributions
   for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
-    list->Add(new TH2F(Form("%s%s", name, fHelper->GetParticleName(idxPart).Data()), 
-                      Form("%s %s Dist %s;Centrality;N_{%s}", title, fHelper->GetParticleTitle(idxPart).Data(), sPtTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), 
-                      24, -0.5, 11.5, 501, -0.5, 500.49));
+    list->Add(new TH3D(Form("h%s%s", name, fHelper->GetParticleName(idxPart).Data()), 
+                      Form("N_{%s} : %s;Centrality;N_{%s}", fHelper->GetParticleTitleLatex(idxPart).Data(), sTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()),
+                      nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 2601, -0.5, 2600.49));
   } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
    
-  list->Add(new TH2F(Form("%sNet%s", name, fHelper->GetParticleName(1).Data()), 
-                    Form("%s Net %s Dist %s;Centrality;N_{%s} - N_{%s}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5, 201, -100.5, 100.49));
-
-  list->Add(new TProfile(Form("%sNet%sM", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s Net %s Dist %s;Centrality;N_{%s} - N_{%s}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()),  24, -0.5, 11.5));
-  list->Add(new TProfile(Form("%sNet%s2M", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s (Net %s)^{2} Dist %s;Centrality;(N_{%s} - N_{%s})^{2}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5));
-  list->Add(new TProfile(Form("%sNet%s3M", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s (Net %s)^{3} Dist %s;Centrality;(N_{%s} - N_{%s})^{3}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5));
-  list->Add(new TProfile(Form("%sNet%s4M", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s (Net %s)^{4} Dist %s;Centrality;(N_{%s} - N_{%s})^{4}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5));
-  list->Add(new TProfile(Form("%sNet%s5M", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s (Net %s)^{5} Dist %s;Centrality;(N_{%s} - N_{%s})^{5}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5));
-  list->Add(new TProfile(Form("%sNet%s6M", name, fHelper->GetParticleName(1).Data()), 
-                        Form("%s (Net %s)^{6} Dist %s;Centrality;(N_{%s} - N_{%s})^{6}", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 24, -0.5, 11.5));
-                        
-  list->Add(new TH2F(Form("%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), 
-                    Form("%s (Net %s)/ Sum Dist %s;Centrality;(N_{%s} - N_{%s})/(N_{%s} + N_{%s})", title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 
-                    24, -0.5, 11.5, 801, -50.5, 50.49));
-
-  if (sName.Contains("fHControlMC")) {
-    for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
-      list->Add(new TH2F(Form("%s%sOver%s", name, fHelper->GetParticleName(idxPart).Data(), fHelper->GetControlParticleName(idxPart).Data()), 
-                        Form("%s %s / %s Dist %s;Centrality;N_{%s}/N_{Tracks}", 
-                             title, fHelper->GetParticleTitle(idxPart).Data(), fHelper->GetControlParticleTitle(idxPart).Data(), sPtTitle.Data(), fHelper->GetParticleTitleLatex(idxPart).Data()), 
-                        24, -0.5, 11.5, 101, 0., 1.));
-    } // for (Int_t idxPart = 0; idxPart < 2; ++idxPart) {
+  // -- Add NetParticle Distributions
+  list->Add(new TH3D(Form("h%sNet%s", name, fHelper->GetParticleName(1).Data()), 
+                    Form("%s : %s;Centrality;%s", sNetTitle.Data(), sTitle.Data(), sNetTitle.Data()),
+                    nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 601, -300.5, 300.49));
+
+  // -- Add NetParticle vs SumParticle
+  list->Add(new TH3D(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()), 
+                    Form("(%s)/(%s) : %s;Centrality;(%s)/(%s)", sNetTitle.Data(), sSumTitle.Data(), sTitle.Data(), sNetTitle.Data(), sSumTitle.Data()), 
+                    nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1], 41, -2.5, 2.49));
+
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <NetParticle^k>
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idx = 1; idx <= fOrder; ++idx) {
+    list->Add(new TProfile2D(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idx), 
+                            Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+                            nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1]));
+  }
+
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <NetParticle^k>  for every SubSample
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) {
+    for (Int_t idx = 1; idx <= fOrder; ++idx) {
+      list->Add(new TProfile2D(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idx, idxSub), 
+                              Form("(%s)^{%d} : %s;Centrality;(%s)^{%d}", sNetTitle.Data(), idx, sTitle.Data(), sNetTitle.Data(), idx),
+                              nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1]));
+    }
   }
+
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <f_ik>
+  // -----------------------------------------------------------------------------------------------
+  list->Add(new TList);
+  TList *fikListPt = static_cast<TList*>(list->Last());
+  fikListPt->SetName(Form("f%sPtFik",name));
+  fikListPt->SetOwner(kTRUE);
+
+  for (Int_t ii = 0; ii <= fOrder; ++ii) {
+    for (Int_t kk = 0; kk <= fOrder; ++kk) {
+      fikListPt->Add(new TProfile2D(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk),
+                                   Form("f_{%02d%02d} : %s;Centrality;#it{p}_{T} (GeV/#it{c});f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), 
+                                   nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1]));
+    }
+  }
+
+  // -- Add counter for number of Non-zero entries
+  for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) 
+    fikListPt->Add(new TH3D(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), idxCent),
+                           Form("f_ik counts : %s Cent %d;i;k;#it{p}_{T} Bins", sTitle.Data(), idxCent),
+                           fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49, nBinsPt+1, -0.5, Double_t(nBinsPt)+0.49));
  
-  //  if (sName.Contains("fHDist")) {
-  //    list->Add(new TH3F(Form("%sNet%s_RecVsMC", name, fHelper->GetParticleName(1).Data()), 
-  //                  Form("%s Net %s Dist - Rec Vs MC - %s;Centrality;(N_{%s} - N_{%s})_{rec};(N_{%s} - N_{%s})_{MC}", 
-  //                       title, fHelper->GetParticleTitle(1).Data(), sPtTitle.Data(),fHelper->GetParticleTitleLatex(1).Data(),fHelper->GetParticleTitleLatex(0).Data()), 
-  //                  24, -0.5, 11.5, 201, -100.5, 100.49, 201, -100.5, 100.49));
-  //  }
-  
+  // -----------------------------------------------------------------------------------------------
+  // -- Add TProfiles for <f_ik> for every SubSample
+  // -----------------------------------------------------------------------------------------------
+  for (Int_t idxSub = 0; idxSub < fHelper->GetNSubSamples(); ++ idxSub) {
+    list->Add(new TList);
+    TList *fikListPtSub = static_cast<TList*>(list->Last());
+    fikListPtSub->SetName(Form("f%sPtFik_%02d",name, idxSub));
+    fikListPtSub->SetOwner(kTRUE);
+    
+    for (Int_t ii = 0; ii <= fOrder; ++ii) {
+      for (Int_t kk = 0; kk <= fOrder; ++kk) {
+       fikListPtSub->Add(new TProfile2D(Form("p%sNet%sF%02d%02d_%02d", name, fHelper->GetParticleName(1).Data(), ii, kk, idxSub),
+                                        Form("f_{%02d%02d} : %s;Centrality;#it{p}_{T} (GeV/#it{c});f_{%02d%02d}", ii, kk, sTitle.Data(), ii, kk), 
+                                        nBinsCent, centBinRange[0], centBinRange[1], nBinsPt, ptBinRange[0], ptBinRange[1]));
+      }
+    }
+
+    // -- Add counter for number of Non-zero entries
+    for (Int_t idxCent = 0; idxCent < nBinsCent; ++idxCent) 
+      fikListPtSub->Add(new TH3D(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), idxCent, idxSub),
+                                Form("f_ik counts : %s Cent %d;i;k;#it{p}_{T} Bins", sTitle.Data(), idxCent),
+                                fOrder+1, -0.5, Double_t(fOrder)+0.49, fOrder+1, -0.5, Double_t(fOrder)+0.49, nBinsPt+1, -0.5, Double_t(nBinsPt)+0.49));
+  }
+
   return;
 }
 
 //________________________________________________________________________
-void AliAnalysisNetParticleDistribution::FillHistSet(const Char_t *name, Float_t *np, Int_t controlIdx)  {
-  // -- Add histogram sets for particle and anti-particle
+void AliAnalysisNetParticleDistribution::FillHistSetCent(const Char_t *name, Int_t idx, Bool_t isMC)  {
+  // -- Fill histogram sets for particle and anti-particle
+  //    dependence : centrality 
+  
+  // -- Get List
+  TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));
+  
+  // -- Get Centrality Bin
+  Float_t centralityBin = fHelper->GetCentralityBin();
+
+  // -- Select MC or Data
+  Int_t **np = (isMC) ? fMCNp : fNp;
 
-  Int_t pid = fHelper->GetParticleSpecies();
-  if( pid < 0 || pid > AliPID::kSPECIES){
-    AliError("Particle ID not in AliPID::kSPECIES");
-    return;
+  // -----------------------------------------------------------------------------------------------
+
+  Int_t sumNp   = np[idx][1]+np[idx][0];  // p + pbar
+  Int_t deltaNp = np[idx][1]-np[idx][0];  // p - pbar
+
+  // -- Fill Particle / Anti-Particle Distributions
+  (static_cast<TH2D*>(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[idx][0]);
+  (static_cast<TH2D*>(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[idx][1]);
+
+  // -- Fill NetParticle Distributions
+  (static_cast<TH2D*>(list->FindObject(Form("h%sNet%s",  name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp);
+
+  // -- Fill NetParticle vs SumParticle
+  Double_t deltaNpOverSumNp = (sumNp == 0.) ? 0. : deltaNp/Double_t(sumNp);
+  (static_cast<TH2D*>(list->FindObject(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpOverSumNp);
+
+  // -----------------------------------------------------------------------------------------------
+
+  Double_t CENT[7] = {1.157990, 1.153218, 1.147382, 1.142218, 1.139568, 1.138152, 1.134517};
+
+  Double_t sumNpX   = np[idx][1]+(np[idx][0]*CENT[Int_t(centralityBin)]);
+  Double_t deltaNpX = np[idx][1]-(np[idx][0]*CENT[Int_t(centralityBin)]);
+
+  // -- Fill Particle / Anti-Particle Distributions
+  (static_cast<TH2D*>(list->FindObject(Form("h%s%sX", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[idx][0]*CENT[Int_t(centralityBin)]);
+  (static_cast<TH2D*>(list->FindObject(Form("h%s%sX", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[idx][1]);
+
+  // -- Fill NetParticle Distributions
+  (static_cast<TH2D*>(list->FindObject(Form("h%sNet%sX",  name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpX);
+
+  // -- Fill NetParticle vs SumParticle
+  Double_t deltaNpXOverSumNpX = (sumNpX == 0.) ? 0. : deltaNpX/sumNpX;
+  (static_cast<TH2D*>(list->FindObject(Form("h%sNet%sOverSumX", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNpXOverSumNpX);
+
+  // -----------------------------------------------------------------------------------------------
+
+  // -- Fill TProfile for <NetParticle^k>
+  Double_t delta = 1.;
+  for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
+    delta *= deltaNp;
+    (static_cast<TProfile*>(list->FindObject(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idxOrder))))->Fill(centralityBin, delta);
+    (static_cast<TProfile*>(list->FindObject(Form("p%sNet%s%dM_%02d", name, fHelper->GetParticleName(1).Data(), idxOrder, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, delta);
   }
 
-  TList *list = static_cast<TList*>(fOutList->FindObject(name));
+  // -- Generate reduced factorials - explictly removing the factorials
+  //    - Reset all to 1
+  for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
+    fRedFactp[idxOrder][0]  = 1.;
+    fRedFactp[idxOrder][1]  = 1.;
+  }
+  
+  //    - start at idx 1, as idx 0 = 1 done by reset
+  for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
+    fRedFactp[idxOrder][0]  = fRedFactp[idxOrder-1][0]  * Double_t(np[idx][0]-(idxOrder-1));
+    fRedFactp[idxOrder][1]  = fRedFactp[idxOrder-1][1]  * Double_t(np[idx][1]-(idxOrder-1));
+  }
 
+  // -- Fill TProfiles for <f_ik> 
+  TList *fikList    = static_cast<TList*>(list->FindObject(Form("f%sFik",name)));
+  TList *fikListSub = static_cast<TList*>(list->FindObject(Form("f%sFik_%02d",name, fHelper->GetSubSampleIdx())));
+  TH2D  *hCntik     = static_cast<TH2D*>(fikList->FindObject(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin))));
+  TH2D  *hCntikSub  = static_cast<TH2D*>(fikListSub->FindObject(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin), fHelper->GetSubSampleIdx())));
+
+  for (Int_t ii = 0; ii <= fOrder; ++ii) {   // ii -> p    -> n1
+    for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
+      // -- use the reduced factorials only 
+      Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];   // n1 *n2 -> p * pbar
+      (static_cast<TProfile*>(fikList->FindObject(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk))))->Fill(centralityBin, fik);
+      (static_cast<TProfile*>(fikListSub->FindObject(Form("p%sNet%sF%02d%02d_%02d", 
+                                                         name, fHelper->GetParticleName(1).Data(), ii, kk, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, fik);
+
+      if (fik != 0.) {
+       hCntik->Fill(ii, kk);
+       hCntikSub->Fill(ii, kk);
+      }
+    }
+  }
+
+  return;
+}
+
+//________________________________________________________________________
+void AliAnalysisNetParticleDistribution::FillHistSetCentPt(const Char_t *name, Int_t idx, Bool_t isMC)  {
+  // -- Add histogram sets for particle and anti-particle
+  //    dependence : centrality and pt
+
+  // -- Get List
+  TList *list = static_cast<TList*>(fOutList->FindObject(Form("f%s",name)));  
+
+  // -- Get Centrality Bin
   Float_t centralityBin = fHelper->GetCentralityBin();
 
-  (static_cast<TH2F*>(list->FindObject(Form("%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[0]);
-  (static_cast<TH2F*>(list->FindObject(Form("%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[1]);
+  // -- Select MC or Data
+  Int_t ***npPt = (isMC) ? fMCNpPt : fNpPt;
+
+  // -----------------------------------------------------------------------------------------------
 
-  Float_t sumNp    = np[1]+np[0];
-  Float_t deltaNp  = np[1]-np[0];
-  Float_t deltaNp2 = deltaNp * deltaNp;
-  Float_t deltaNp3 = deltaNp2 * deltaNp;
+  // -- Loop over the pt bins
+  for (Int_t idxPt  = 0; idxPt < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++idxPt) {
+    
+    Int_t deltaNp = npPt[idx][1][idxPt]-npPt[idx][0][idxPt]; // p - pbar
+    Int_t sumNp   = npPt[idx][1][idxPt]+npPt[idx][0][idxPt]; // p + pbar
 
-  (static_cast<TH2F*>(list->FindObject(Form("%sNet%s",  name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp);
+    // -- Fill Particle / Anti-Particle Distributions
+    (static_cast<TH3D*>(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, idxPt, npPt[idx][0][idxPt]);
+    (static_cast<TH3D*>(list->FindObject(Form("h%s%s", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, npPt[idx][1][idxPt]);
+    
+    // -- Fill NetParticle Distributions
+    (static_cast<TH3D*>(list->FindObject(Form("h%sNet%s",  name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, deltaNp);
+    
+    // -- Fill NetParticle vs SumParticle
+    Double_t deltaNpOverSumNp = (sumNp == 0.) ? 0. : deltaNp/Double_t(sumNp);
+    (static_cast<TH3D*>(list->FindObject(Form("h%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, idxPt, deltaNpOverSumNp);
 
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%sM",  name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp);
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%s2M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2);
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%s3M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2*deltaNp);
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%s4M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp2*deltaNp2);
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%s5M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp3*deltaNp2);
-  (static_cast<TProfile*>(list->FindObject(Form("%sNet%s6M", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp3*deltaNp3);
+    // -----------------------------------------------------------------------------------------------
 
-  (static_cast<TH2F*>(list->FindObject(Form("%sNet%sOverSum", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, deltaNp/sumNp);
+    // -- Fill TProfile for <NetParticle^k>
+    Double_t delta = 1.;
+    for (Int_t idxOrder = 1; idxOrder <= fOrder; ++idxOrder) {
+      delta *= deltaNp;
+      (static_cast<TProfile2D*>(list->FindObject(Form("p%sNet%s%dM", name, fHelper->GetParticleName(1).Data(), idxOrder))))->Fill(centralityBin, idxPt, delta);
+      (static_cast<TProfile2D*>(list->FindObject(Form("p%sNet%s%dM_%02d", 
+                                                     name, fHelper->GetParticleName(1).Data(), idxOrder, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, idxPt, delta);
+    }
+    
+    // -- Generate reduced factorials - explictly removing the factorials
+    //    - Reset all to 1
+    for (Int_t idxOrder = 0; idxOrder <= fOrder; ++ idxOrder) {
+      fRedFactp[idxOrder][0] = 1.;
+      fRedFactp[idxOrder][1] = 1.;
+    }
 
-  TString sName(name);
-  if (sName.Contains("fHControlMC") && controlIdx >= 0) {
-    (static_cast<TH2F*>(list->FindObject(Form("%s%sOverlambdabar", name, fHelper->GetParticleName(0).Data()))))->Fill(centralityBin, np[0]/fControlMCNp[controlIdx][0]);
-    (static_cast<TH2F*>(list->FindObject(Form("%s%sOverlambda", name, fHelper->GetParticleName(1).Data()))))->Fill(centralityBin, np[1]/fControlMCNp[controlIdx][1]);
-  }
+    //    - start at idx 1, as idx 0 = 1 done by reset
+    for (Int_t idxOrder = 1; idxOrder <= fOrder; ++ idxOrder) {
+      fRedFactp[idxOrder][0] = fRedFactp[idxOrder-1][0] * Double_t(npPt[idx][0][idxPt]-(idxOrder-1));
+      fRedFactp[idxOrder][1] = fRedFactp[idxOrder-1][1] * Double_t(npPt[idx][1][idxPt]-(idxOrder-1));
+    }
 
-  //  if (sName.Contains("fHDist")) {
-  //    Float_t deltaMCNp  = fMCNp[controlIdx][1]-fMCNp[controldx][0];
-  //    (static_cast<TH3F*>(list->FindObject(Form("%sNetp_RecVsMC", name))))->Fill(centralityBin, deltaNp, deltaMCNp);
-  //  }
+    // -- Fill TProfiles for <f_ik> 
+    TList *fikListPt    = static_cast<TList*>(list->FindObject(Form("f%sPtFik",name)));
+    TList *fikListPtSub = static_cast<TList*>(list->FindObject(Form("f%sPtFik_%02d",name, fHelper->GetSubSampleIdx())));
+    TH3D  *hCntikPt     = static_cast<TH3D*>(fikListPt->FindObject(Form("p%sNet%sFCounts_%02d", name, fHelper->GetParticleName(1).Data(), Int_t(centralityBin))));
+    TH3D  *hCntikPtSub  = static_cast<TH3D*>(fikListPtSub->FindObject(Form("p%sNet%sFCounts_%02d_%02d", name, fHelper->GetParticleName(1).Data(), 
+                                                                          Int_t(centralityBin), fHelper->GetSubSampleIdx())));
+    for (Int_t ii = 0; ii <= fOrder; ++ii) {   // ii -> p    -> n1
+      for (Int_t kk = 0; kk <= fOrder; ++kk) { // kk -> pbar -> n2
+       Double_t fik = fRedFactp[ii][1] * fRedFactp[kk][0];   // n1 *n2 -> p * pbar
+       (static_cast<TProfile2D*>(fikListPt->FindObject(Form("p%sNet%sF%02d%02d", name, fHelper->GetParticleName(1).Data(), ii, kk))))->Fill(centralityBin, idxPt, fik);
+       (static_cast<TProfile2D*>(fikListPtSub->FindObject(Form("p%sNet%sF%02d%02d_%02d", 
+                                                               name, fHelper->GetParticleName(1).Data(), ii, kk, fHelper->GetSubSampleIdx()))))->Fill(centralityBin, idxPt, fik);
+       
+       if (fik != 0.) {
+         hCntikPt->Fill(ii, kk, idxPt);
+         hCntikPtSub->Fill(ii, kk, idxPt);
+       }
+      }
+    }
+    
+  } // for (Int_t idxPt  = 0; idxPt < AliAnalysisNetParticleHelper::fgkfHistNBinsPt; ++idxPt) {
 
   return;
 }