diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h
index 45e4f7c26d7..9623f795ea6 100644
--- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h
+++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h
@@ -98,6 +98,7 @@ struct Pi0EtaToGammaGamma {
   o2::framework::ConfigurableAxis ConfCentBins{"ConfCentBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"};
   o2::framework::ConfigurableAxis ConfEPBins{"ConfEPBins", {o2::framework::VARIABLE_WIDTH, -o2::constants::math::PIHalf, -o2::constants::math::PIQuarter, 0.0f, +o2::constants::math::PIQuarter, +o2::constants::math::PIHalf}, "Mixing bins - event plane angle"};
   o2::framework::ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {o2::framework::VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"};
+  o2::framework::Configurable<bool> cfgSparseFullAxes{"cfgSparseFullAxes", false, "add event plane, occupancy and z-vertex axes to Pair/same/hs and Pair/mix/hs"};
 
   o2::framework::Configurable<int> cfgAlphaMesonCut{"cfgAlphaMesonCut", 0, "flag for photon energy asymmetry distribution cut: 0: no cut, 1: cut specific value, 2: cut depending on pT"};
   o2::framework::Configurable<float> cfgAlphaMeson{"cfgAlphaMeson", 0.65, "photon energy asymmetry distribution parameter for specific value cut"};
@@ -423,10 +424,13 @@ struct Pi0EtaToGammaGamma {
     emh2 = new o2::aod::pwgem::dilepton::utils::EventMixingHandler<std::tuple<int, int, int, int>, std::pair<int, int>, o2::aod::pwgem::photonmeson::utils::EMPhoton>(ndepth);
 
     o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(&fRegistry);
+    std::vector<double> occBinsForHs(occ_bin_edges.begin(), occ_bin_edges.end());
+    std::vector<double> epBinsForHs(ep_bin_edges.begin(), ep_bin_edges.end());
+    std::vector<double> vtxBinsForHs(zvtx_bin_edges.begin(), zvtx_bin_edges.end());
     if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) {
-      o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, "ee#gamma");
+      o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs, "ee#gamma");
     } else {
-      o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, "#gamma#gamma");
+      o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs, "#gamma#gamma");
     }
     DefineEMEventCut();
     DefinePCMCut();
@@ -800,6 +804,8 @@ struct Pi0EtaToGammaGamma {
         occbin = static_cast<int>(occ_bin_edges.size()) - 2;
       }
 
+      float occupancy = (cfgOccupancyEstimator == 1) ? static_cast<float>(collision.trackOccupancyInTimeRange()) : collision.ft0cOccupancyInTimeRange();
+
       // LOGF(info, "collision.globalIndex() = %d, collision.posZ() = %f, centrality = %f, ep2 = %f, collision.trackOccupancyInTimeRange() = %d, zbin = %d, centbin = %d, epbin = %d, occbin = %d", collision.globalIndex(), collision.posZ(), centrality, ep2, collision.trackOccupancyInTimeRange(), zbin, centbin, epbin, occbin);
 
       std::tuple<int, int, int, int> key_bin = std::make_tuple(zbin, centbin, epbin, occbin);
@@ -859,7 +865,11 @@ struct Pi0EtaToGammaGamma {
               continue;
             }
 
-            fRegistry.fill(HIST("Pair/same/hs"), veeg.M(), veeg.Pt(), weight);
+            if (cfgSparseFullAxes.value) {
+              fRegistry.fill(HIST("Pair/same/hs"), veeg.M(), veeg.Pt(), ep2, occupancy, collision.posZ(), weight);
+            } else {
+              fRegistry.fill(HIST("Pair/same/hs"), veeg.M(), veeg.Pt(), weight);
+            }
 
             std::pair<int, int> tuple_tmp_id2 = std::make_pair(pos2.trackId(), ele2.trackId());
             if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g1.globalIndex()) == used_photonIds_per_col.end()) {
@@ -953,7 +963,11 @@ struct Pi0EtaToGammaGamma {
             wpair *= g2.omegaMBWeight();
           }
 
-          fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), wpair);
+          if (cfgSparseFullAxes.value) {
+            fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), ep2, occupancy, collision.posZ(), wpair);
+          } else {
+            fRegistry.fill(HIST("Pair/same/hs"), v12.M(), v12.Pt(), wpair);
+          }
 
           if (std::find(used_photonIds_per_col.begin(), used_photonIds_per_col.end(), g1.globalIndex()) == used_photonIds_per_col.end()) {
             auto emphoton1 = o2::aod::pwgem::photonmeson::utils::EMPhoton(g1.pt(), g1.eta(), g1.phi(), 0);
@@ -1047,7 +1061,11 @@ struct Pi0EtaToGammaGamma {
                 continue;
               }
 
-              fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              if (cfgSparseFullAxes.value) {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), ep2, occupancy, collision.posZ(), weight);
+              } else {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              }
             }
           }
         } // end of loop over mixed event pool
@@ -1082,7 +1100,11 @@ struct Pi0EtaToGammaGamma {
               if (std::fabs(v12.Rapidity()) > maxY) {
                 continue;
               }
-              fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              if (cfgSparseFullAxes.value) {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), ep2, occupancy, collision.posZ(), weight);
+              } else {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              }
             }
           }
         } // end of loop over mixed event pool
@@ -1115,7 +1137,11 @@ struct Pi0EtaToGammaGamma {
               if (std::fabs(v12.Rapidity()) > maxY) {
                 continue;
               }
-              fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              if (cfgSparseFullAxes.value) {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), ep2, occupancy, collision.posZ(), weight);
+              } else {
+                fRegistry.fill(HIST("Pair/mix/hs"), v12.M(), v12.Pt(), weight);
+              }
             }
           }
         } // end of loop over mixed event pool
diff --git a/PWGEM/PhotonMeson/Utils/NMHistograms.h b/PWGEM/PhotonMeson/Utils/NMHistograms.h
index ac361d9783e..81a5d712669 100644
--- a/PWGEM/PhotonMeson/Utils/NMHistograms.h
+++ b/PWGEM/PhotonMeson/Utils/NMHistograms.h
@@ -32,7 +32,12 @@
 
 namespace o2::aod::pwgem::photonmeson::utils::nmhistogram
 {
-inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool isMC, const char* pairname = "#gamma#gamma")
+inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool isMC, 
+  bool sparseFullAxes = false,
+  std::vector<double> const& occBins = {0, 500, 1000, 2000, 3000, 5000, 10000},
+  std::vector<double> const& epBins  = {-M_PI/2, -M_PI/4, 0, M_PI/4, M_PI/2},
+  std::vector<double> const& vtxBins = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10},
+  const char* pairname = "#gamma#gamma")
 {
   // !!Don't change pt,eta,y binning. These binnings have to be consistent with binned data at skimming.!!
   std::vector<double> ptbins;
@@ -51,6 +56,9 @@ inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool is
   }
   const o2::framework::AxisSpec axis_pt{ptbins, Form("p_{T,%s} (GeV/c)", pairname)};
   const o2::framework::AxisSpec axis_mass{400, 0, 0.8, Form("m_{%s} (GeV/c^{2})", pairname)};
+  const o2::framework::AxisSpec axis_occ{occBins, "occupancy"};
+  const o2::framework::AxisSpec axis_ep {epBins,  "#Psi_{2} (rad)"};
+  const o2::framework::AxisSpec axis_vtx{vtxBins, "Z_{vtx} (cm)"};
 
   if (isMC) {
     fRegistry->add("Pair/Pi0/hs_Primary", "rec. true pi0", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt}, true);
@@ -74,9 +82,19 @@ inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool is
     fRegistry->get<TH2>(HIST("Generated/Eta/hPtY"))->SetXTitle("p_{T} (GeV/c)");
     fRegistry->get<TH2>(HIST("Generated/Eta/hPtY"))->SetYTitle("rapidity |y|");
   } else {
-    fRegistry->add("Pair/same/hs", "diphoton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt}, true);
+    if (sparseFullAxes) {
+      fRegistry->add("Pair/same/hs", "diphoton",
+            o2::framework::HistType::kTHnSparseD,
+            {axis_mass, axis_pt, axis_ep, axis_occ, axis_vtx}, true);
+    } else {
+      fRegistry->add("Pair/same/hs", "diphoton",
+            o2::framework::HistType::kTHnSparseD,
+            {axis_mass, axis_pt}, true);
+    }
     fRegistry->addClone("Pair/same/", "Pair/mix/");
   }
+
+  
 }
 
 template <typename TDiphoton, o2::soa::is_iterator TMCParitlce, o2::soa::is_table TMCParticles, o2::soa::is_table TMCCollisions>