diff --git a/PhysicsTools/NanoAOD/plugins/NPUTablesProducer.cc b/PhysicsTools/NanoAOD/plugins/NPUTablesProducer.cc
index d8bd6999647df..aebb3ac83b200 100644
--- a/PhysicsTools/NanoAOD/plugins/NPUTablesProducer.cc
+++ b/PhysicsTools/NanoAOD/plugins/NPUTablesProducer.cc
@@ -10,13 +10,15 @@
 
 #include <vector>
 #include <iostream>
+#include <algorithm>
 
 class NPUTablesProducer : public edm::global::EDProducer<> {
 public:
   NPUTablesProducer(edm::ParameterSet const& params)
       : npuTag_(consumes<std::vector<PileupSummaryInfo>>(params.getParameter<edm::InputTag>("src"))),
         pvTag_(consumes<std::vector<reco::Vertex>>(params.getParameter<edm::InputTag>("pvsrc"))),
-        vz_(params.getParameter<std::vector<double>>("zbins")) {
+        vz_(params.getParameter<std::vector<double>>("zbins")),
+        savePtHatMax_(params.getParameter<bool>("savePtHatMax")) {
     produces<nanoaod::FlatTable>();
   }
 
@@ -47,6 +49,8 @@ class NPUTablesProducer : public edm::global::EDProducer<> {
     float pudensity = 0;
     float gpudensity = 0;
 
+    float pthatmax = 0;
+
     for (unsigned int ibx = 0; ibx < npuProd.size(); ibx++) {
       if (npuProd[ibx].getBunchCrossing() == 0) {
         bx0 = ibx;
@@ -64,6 +68,10 @@ class NPUTablesProducer : public edm::global::EDProducer<> {
             gpudensity++;
         }
         gpudensity /= (20.0 * (*(zbin) - *(zbin - 1)));
+
+        if (savePtHatMax_) {
+          pthatmax = *max_element(npuProd[ibx].getPU_pT_hats().begin(), npuProd[ibx].getPU_pT_hats().end());
+        }
       }
     }
     unsigned int eoot = 0;
@@ -86,6 +94,9 @@ class NPUTablesProducer : public edm::global::EDProducer<> {
     out.addColumnValue<int>("sumLOOT", loot, "number of late out of time pileup");
     out.addColumnValue<float>("pudensity", pudensity, "PU vertices / mm");
     out.addColumnValue<float>("gpudensity", gpudensity, "Generator-level PU vertices / mm");
+    if (savePtHatMax_) {
+      out.addColumnValue<float>("pthatmax", pthatmax, "Maximum pt-hat");
+    }
   }
 
   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
@@ -95,6 +106,7 @@ class NPUTablesProducer : public edm::global::EDProducer<> {
     desc.add<edm::InputTag>("pvsrc", edm::InputTag("offlineSlimmedPrimaryVertices"))->setComment("tag for the PVs");
     desc.add<std::vector<double>>("zbins", {})
         ->setComment("Z bins to compute the generator-level number of PU vertices per mm");
+    desc.add<bool>("savePtHatMax", false)->setComment("Store maximum pt-hat of PU");
     descriptions.add("puTable", desc);
   }
 
@@ -103,6 +115,8 @@ class NPUTablesProducer : public edm::global::EDProducer<> {
   const edm::EDGetTokenT<std::vector<reco::Vertex>> pvTag_;
 
   const std::vector<double> vz_;
+
+  bool savePtHatMax_;
 };
 
 #include "FWCore/Framework/interface/MakerMacros.h"
diff --git a/PhysicsTools/NanoAOD/python/custom_jme_cff.py b/PhysicsTools/NanoAOD/python/custom_jme_cff.py
index 1829e4c538466..38222859ede37 100644
--- a/PhysicsTools/NanoAOD/python/custom_jme_cff.py
+++ b/PhysicsTools/NanoAOD/python/custom_jme_cff.py
@@ -874,6 +874,12 @@ def PrepJMECustomNanoAOD(process,runOnMC):
     cfg = { k : v for k, v in jetConfig.items() if k != "enabled"}
     recoJetInfo = recoJA.addRecoJetCollection(process, **cfg)
     AddNewPatJets(process, recoJetInfo, runOnMC)
+
+  ###########################################################################
+  # Save Maximum of Pt Hat Max
+  ###########################################################################
+  if runOnMC:
+    process.puTable.savePtHatMax = True
   
   return process
 
diff --git a/PhysicsTools/NanoAOD/python/globals_cff.py b/PhysicsTools/NanoAOD/python/globals_cff.py
index bfa17e721636b..acaf0711e4175 100644
--- a/PhysicsTools/NanoAOD/python/globals_cff.py
+++ b/PhysicsTools/NanoAOD/python/globals_cff.py
@@ -14,7 +14,8 @@
 puTable = cms.EDProducer("NPUTablesProducer",
         src = cms.InputTag("slimmedAddPileupInfo"),
         pvsrc = cms.InputTag("offlineSlimmedPrimaryVertices"),
-        zbins = cms.vdouble( [0.0,1.7,2.6,3.0,3.5,4.2,5.2,6.0,7.5,9.0,12.0] )
+        zbins = cms.vdouble( [0.0,1.7,2.6,3.0,3.5,4.2,5.2,6.0,7.5,9.0,12.0] ),
+        savePtHatMax = cms.bool(False), 
 )
 
 genTable  = cms.EDProducer("SimpleGenEventFlatTableProducer",