From 92c3aa22fbb6ae2309572ce69920ade5f3f53124 Mon Sep 17 00:00:00 2001
From: kevScheuer <kscheuer@jlab.org>
Date: Wed, 30 Sep 2026 10:15:38 -0400
Subject: [PATCH 1/2] Implement poisson bootstrapping

---
 AmpTools/GPUManager/GPUManager.cc             |  9 +++-
 AmpTools/GPUManager/GPUManager.h              |  1 +
 AmpTools/IUAmpTools/AmpToolsInterface.cc      | 38 +++++++++++++++-
 AmpTools/IUAmpTools/AmpToolsInterface.h       |  9 ++++
 AmpTools/IUAmpTools/AmpVecs.cc                | 36 ++++++++++++++--
 AmpTools/IUAmpTools/AmpVecs.h                 | 22 +++++++++-
 AmpTools/IUAmpTools/AmplitudeManager.cc       |  6 ++-
 AmpTools/IUAmpTools/LikelihoodCalculator.cc   | 24 +++++++++++
 AmpTools/IUAmpTools/LikelihoodCalculator.h    |  4 +-
 AmpTools/IUAmpTools/NormIntInterface.cc       | 16 +++++++
 AmpTools/IUAmpTools/NormIntInterface.h        |  5 ++-
 .../IUAmpToolsMPI/LikelihoodCalculatorMPI.cc  | 29 +++++++++++++
 .../IUAmpToolsMPI/LikelihoodCalculatorMPI.h   |  4 ++
 .../IUAmpToolsMPI/LikelihoodManagerMPI.cc     |  5 +++
 AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h |  1 +
 AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc | 43 ++++++++++++-------
 AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.h  |  2 +
 17 files changed, 225 insertions(+), 29 deletions(-)

diff --git a/AmpTools/GPUManager/GPUManager.cc b/AmpTools/GPUManager/GPUManager.cc
index 05dc7009..92c42713 100644
--- a/AmpTools/GPUManager/GPUManager.cc
+++ b/AmpTools/GPUManager/GPUManager.cc
@@ -345,8 +345,7 @@ GPUManager::copyDataToGPU( const AmpVecs& a, bool use4Vectors )
   #endif
   
   // copy the weights to the GPU
-  gpuErrChk( cudaMemcpy( m_pfDevWeights, a.m_pdWeights,
-                         m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) );
+  copyWeightsToGPU( a );
   
   // we only need to copy the four vectors to the GPU if the
   // amplitude evaulation kernels need them -- this is done by
@@ -426,6 +425,12 @@ GPUManager::copyUserVarsToGPU( const AmpVecs& a )
   #endif
 }
 
+void
+GPUManager::copyWeightsToGPU( const AmpVecs& a )
+{
+  gpuErrChk( cudaMemcpy( m_pfDevWeights, a.m_pdWeights,
+                         m_iGDoubleDataArrSize, cudaMemcpyHostToDevice ) );
+}
 
 void
 GPUManager::copyAmpsFromGPU( AmpVecs& a )
diff --git a/AmpTools/GPUManager/GPUManager.h b/AmpTools/GPUManager/GPUManager.h
index c5d64934..056cc3c7 100644
--- a/AmpTools/GPUManager/GPUManager.h
+++ b/AmpTools/GPUManager/GPUManager.h
@@ -80,6 +80,7 @@ class GPUManager
   
   void copyDataToGPU( const AmpVecs& a, bool use4Vectors = true  );
   void copyUserVarsToGPU( const AmpVecs& a );
+  void copyWeightsToGPU( const AmpVecs& a );
   
   void calcAmplitudeAll( const Amplitude* amp, size_t uAmpFactOffset,
                          const vector< vector< int > >* pvPermutations,
diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.cc b/AmpTools/IUAmpTools/AmpToolsInterface.cc
index d4cb3c49..2561c54d 100644
--- a/AmpTools/IUAmpTools/AmpToolsInterface.cc
+++ b/AmpTools/IUAmpTools/AmpToolsInterface.cc
@@ -33,7 +33,8 @@
 // held liable for any liability with respect to any claim by the user or
 // any other party arising from use of the program.
 //******************************************************************************
-
+#include <ctime>
+#include <map>
 
 #include "IUAmpTools/AmpToolsInterface.h"
 #include "MinuitInterface/MinuitMinimizationManager.h"
@@ -400,6 +401,41 @@ AmpToolsInterface::randomizeParameter( const string& parName, float min, float m
   invalidateAmps();
 }
 
+void
+AmpToolsInterface::bootstrap( unsigned int seed ){
+  
+  if( m_functionality != kFull ) return;
+  if ( seed == 0 ) seed = (unsigned int) time( NULL );
+
+  // give shared acc MC datasets across reactions the same bootstrap seed. Each will
+  // still draw and reassign its own set of weights, but they will be identical. Gen MC 
+  // datasets are never bootstrapped.
+  map< DataReader*, unsigned int > mcSeedForReader;
+  unsigned int nGroups = 0;
+
+  for (unsigned int irct = 0; irct < m_configurationInfo->reactionList().size(); irct++){
+    
+    ReactionInfo* reaction = m_configurationInfo->reactionList()[irct];
+    string reactionName( reaction->reactionName() );
+
+    unsigned int mcSeed = 0;
+    if( DataReader* accMC = accMCReader( reactionName ) ){
+      map< DataReader*, unsigned int >::iterator it = mcSeedForReader.find( accMC );
+      if( it == mcSeedForReader.end() ){
+        mcSeed = seed + 1000000u * (++nGroups);
+        mcSeedForReader[accMC] = mcSeed;
+      }
+      else mcSeed = it->second;
+    }
+
+    // data/bkgnd files do not share common seeds across reactions
+    unsigned int dataSeed = seed + irct;
+
+    if (LikelihoodCalculator* likCalc = likelihoodCalculator( reactionName ) )
+      likCalc->applyPoissonWeights( dataSeed, mcSeed ); 
+  }
+}
+
 void
 AmpToolsInterface::finalizeFit( const string& tag ){
   
diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.h b/AmpTools/IUAmpTools/AmpToolsInterface.h
index 24c7b8af..ac7dc6e2 100644
--- a/AmpTools/IUAmpTools/AmpToolsInterface.h
+++ b/AmpTools/IUAmpTools/AmpToolsInterface.h
@@ -259,6 +259,15 @@ class AmpToolsInterface{
    */
   
   void randomizeParameter( const string& parName, float min = 0, float max = 1 );
+
+  /** This function will perform Poisson bootstrapping by re-weighting the 
+   * signal, accMC, and background (if available) events, using the provided 
+   * seed. Data and background are always given the same seed per reaction, while
+   * accMC has its own unique, but still traceable, seed. Repeated accMC samples
+   * use the same seed.
+   */
+
+  void bootstrap( unsigned int seed = 0 );
   
   /** Print final fit results to a file.  The tag can be used to
    *  generate a unique name in the case that multiple results are
diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc
index 4ea5c40a..da218428 100644
--- a/AmpTools/IUAmpTools/AmpVecs.cc
+++ b/AmpTools/IUAmpTools/AmpVecs.cc
@@ -59,8 +59,9 @@ AmpVecs::AmpVecs(){
   m_maxFactPerEvent = 0 ;
   m_userVarsPerEvent = 0;
   
-  m_pdData      = 0 ;
-  m_pdWeights   = 0 ;
+  m_pdData            = 0 ;
+  m_pdWeights         = 0 ;
+  m_pdOriginalWeights = 0 ;
   
   m_pdAmps       = 0 ;
   m_pdAmpFactors = 0 ;
@@ -104,6 +105,10 @@ AmpVecs::deallocAmpVecs()
   if(m_pdWeights)
     delete[] m_pdWeights;
   m_pdWeights=0;
+
+  if(m_pdOriginalWeights)
+    delete[] m_pdOriginalWeights;
+  m_pdOriginalWeights=0;
   
   deallocTerms( true );
 
@@ -206,6 +211,7 @@ AmpVecs::loadEvent( const Kinematics* pKinematics, size_t iEvent,
 
     m_pdData = new GDouble[4*m_iNParticles*m_iNEvents];
     m_pdWeights = new GDouble[m_iNEvents];        
+    m_pdOriginalWeights = new GDouble[m_iNEvents];
   }
     
   loadDataArrayElement( pKinematics, iEvent );
@@ -222,7 +228,7 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly, size_t chunk
   
   //  Make sure no data is already loaded
   
-  if( m_pdData != NULL || m_pdWeights != NULL ){
+  if( m_pdData != NULL || m_pdWeights != NULL || m_pdOriginalWeights != NULL ){
 
     report( ERROR, kModule ) << "Trying to load data into a non-empty AmpVecs object\n"<<flush;
     assert(false);
@@ -276,6 +282,7 @@ AmpVecs::loadData( DataReader* pDataReader, bool needsUserVarsOnly, size_t chunk
       m_iNParticles = pKinematics->particleList().size();
       m_pdData = new GDouble[4*m_iNParticles*m_iNEvents];
       m_pdWeights = new GDouble[m_iNEvents]; 
+      m_pdOriginalWeights = new GDouble[m_iNEvents];
 
       #ifdef GPU_ACCELERATION
 
@@ -337,6 +344,26 @@ AmpVecs::loadDataArrayElement( const Kinematics* pKinematics, size_t iEvent ){
   }
 
   m_pdWeights[iEvent] = pKinematics->weight();
+  m_pdOriginalWeights[iEvent] = pKinematics->weight();
+}
+
+void
+AmpVecs::applyPoissonWeights( const int* weights ){
+
+  m_dSumWeights = 0;
+
+  for(size_t iEvent = 0; iEvent < m_iNTrueEvents; iEvent++){
+
+    m_pdWeights[iEvent] = m_pdOriginalWeights[iEvent] * weights[iEvent];
+    m_dSumWeights += m_pdWeights[iEvent];
+  }
+
+  m_integralValid = false;
+
+  #ifdef GPU_ACCELERATION
+    m_gpuMan.copyWeightsToGPU( *this );
+  #endif
+
 }
 
 void
@@ -522,6 +549,9 @@ AmpVecs::shareDataWith( AmpVecs* targetAmpVecs, bool needsUserVarsOnly ){
   targetAmpVecs->m_pdWeights = new GDouble[m_iNEvents];
   memcpy( targetAmpVecs->m_pdWeights, m_pdWeights, 
 	  sizeof(GDouble)*m_iNEvents );
+  targetAmpVecs->m_pdOriginalWeights = new GDouble[m_iNEvents];
+  memcpy( targetAmpVecs->m_pdOriginalWeights, m_pdOriginalWeights,
+    sizeof(GDouble)*m_iNEvents );
   
   targetAmpVecs->m_dataLoaded = true;
   targetAmpVecs->m_hasNonUnityWeights = m_hasNonUnityWeights;
diff --git a/AmpTools/IUAmpTools/AmpVecs.h b/AmpTools/IUAmpTools/AmpVecs.h
index d78a5cd3..0916b134 100644
--- a/AmpTools/IUAmpTools/AmpVecs.h
+++ b/AmpTools/IUAmpTools/AmpVecs.h
@@ -81,7 +81,7 @@ struct AmpVecs
   size_t m_iNTrueEvents;
   
   /**
-   * A double that stores the absolute value of the sum of the weights.  (For
+   * A double that stores the absolute value of the sum of the current weights.  (For
    * cases where all weights are unity, this is simply the number of true 
    * events.)
    */
@@ -118,10 +118,16 @@ struct AmpVecs
   GDouble* m_pdData;
   
   /**
-   * An array of length iNEvents that stores the event weights for each event.
+
+  * An array of length iNEvents that stores the current event weights for each event.
    */
   GDouble* m_pdWeights;
   
+  /**
+   * An array of length iNEvents that stores the original event weights for each event.
+   */
+  GDouble* m_pdOriginalWeights;
+  
   /**
    * An array of length 2 * iNAmps * iNEvents that stores the real and imaginary
    * parts of the complete decay amplitude (product of factors) for each event.
@@ -296,6 +302,18 @@ struct AmpVecs
    */
   void loadEvent( const Kinematics* pKinematics, size_t iEvent = 0,
                   size_t iNTrueEvents = 1, bool needsUserVarsOnly = false );
+
+  /**
+   * This will overwrite the m_pdWeights array with the original weights multiplied by 
+   * an integer randomly drawn from a Poisson distribution of mean 1. For any given 
+   * weighted event, this will effectively
+   * - (integer == 0) remove the event 
+   * - (integer == 1) keep the event 
+   * - (integer == N) repeat the event N times, for N > 1
+   * 
+   * \param[in] weights an array of integers of length m_iNTrueEvents
+   */
+  void applyPoissonWeights(const int* weights);
   
   /**
    * A helper routine to get an event i from the array of data and weights.
diff --git a/AmpTools/IUAmpTools/AmplitudeManager.cc b/AmpTools/IUAmpTools/AmplitudeManager.cc
index c5522235..104b9d69 100644
--- a/AmpTools/IUAmpTools/AmplitudeManager.cc
+++ b/AmpTools/IUAmpTools/AmplitudeManager.cc
@@ -744,14 +744,16 @@ SCOREP_USER_REGION_BEGIN( calcSumLogIntensity, "calcSumLogIntensity", SCOREP_USE
   calcIntensities( a );
   
   for( int iEvent=0; iEvent < a.m_iNTrueEvents; iEvent++ ){
+
+    double weight = a.m_pdWeights[iEvent];
+    if( weight == 0 ) continue; // avoid divide by zero
     
     // here divide out the weight that was put into the intensity calculation
     // and weight the log -- in practice this just contributes an extra constant
     // term in the likelihood equal to sum -w_i * log( w_i ), but the division
     // helps avoid problems with negative weights, which may be used
     // in background subtraction
-    dSumLogI += a.m_pdWeights[iEvent] *
-    G_LOG( a.m_pdIntensity[iEvent] / a.m_pdWeights[iEvent] );
+    dSumLogI += weight * G_LOG( a.m_pdIntensity[iEvent] / weight );
   }
   
 #else
diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.cc b/AmpTools/IUAmpTools/LikelihoodCalculator.cc
index af659398..f0ebb74a 100644
--- a/AmpTools/IUAmpTools/LikelihoodCalculator.cc
+++ b/AmpTools/IUAmpTools/LikelihoodCalculator.cc
@@ -38,6 +38,8 @@
 #include <cassert>
 #include <string>
 
+#include "TRandom3.h"
+
 #include "IUAmpTools/LikelihoodCalculator.h"
 #include "IUAmpTools/IntensityManager.h"
 #include "IUAmpTools/DataReader.h"
@@ -305,6 +307,28 @@ SCOREP_USER_REGION_END( dataTerm )
   return sumLnI;
 }
 
+void
+LikelihoodCalculator::applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed ){
+  if( m_firstDataCalc ) dataTerm();
+
+  TRandom3 rndGen( dataSeed );
+
+  vector<int> sigW( m_ampVecsSignal.m_iNTrueEvents);
+  for( size_t i = 0; i < sigW.size(); ++i) sigW[i] = rndGen.Poisson( 1.0 );
+  m_ampVecsSignal.applyPoissonWeights( &sigW[0] );
+  m_sumDataWeights = m_ampVecsSignal.m_dSumWeights;
+
+  if( m_hasBackground ){
+    vector<int> bkgW( m_ampVecsBkgnd.m_iNTrueEvents);
+    for( size_t i = 0; i < bkgW.size(); ++i) bkgW[i] = rndGen.Poisson( 1.0 );
+    m_ampVecsBkgnd.applyPoissonWeights( &bkgW[0] );
+    m_sumBkgWeights = m_ampVecsBkgnd.m_dSumWeights;
+  }
+
+  m_normInt.applyPoissonWeights( mcSeed );
+  
+}
+
 void
 LikelihoodCalculator::invalidateTerms(){
   
diff --git a/AmpTools/IUAmpTools/LikelihoodCalculator.h b/AmpTools/IUAmpTools/LikelihoodCalculator.h
index 62df8c0d..cb59a68a 100644
--- a/AmpTools/IUAmpTools/LikelihoodCalculator.h
+++ b/AmpTools/IUAmpTools/LikelihoodCalculator.h
@@ -85,6 +85,7 @@ class LikelihoodCalculator : public MIFunctionContribution
   double operator()();
   
   virtual double numSignalEvents();
+  virtual void applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed );
   
   void invalidateTerms();
   
@@ -108,13 +109,14 @@ class LikelihoodCalculator : public MIFunctionContribution
   void setNumBkgEvents (  double num ) { m_numBkgEvents   = num; }
   void setSumDataWeights(  double sum ) { m_sumDataWeights = sum; }
   void setNumDataEvents(  double num ) { m_numDataEvents  = num; }
+
+  const NormIntInterface& m_normInt;
   
 private:
   
   bool m_hasBackground;
   
   const IntensityManager& m_intenManager;
-  const NormIntInterface& m_normInt;
 
   DataReader* m_dataReaderSignal;
   DataReader* m_dataReaderBkgnd;
diff --git a/AmpTools/IUAmpTools/NormIntInterface.cc b/AmpTools/IUAmpTools/NormIntInterface.cc
index f4fd296b..dacf904c 100644
--- a/AmpTools/IUAmpTools/NormIntInterface.cc
+++ b/AmpTools/IUAmpTools/NormIntInterface.cc
@@ -44,6 +44,8 @@
 #include <cstring>
 #include <iomanip>
 
+#include "TRandom3.h"
+
 #include "IUAmpTools/NormIntInterface.h"
 #include "IUAmpTools/report.h"
 const char* NormIntInterface::kModule = "NormIntInterface";
@@ -603,6 +605,20 @@ NormIntInterface::loadMC() const {
   m_sumGenWeights = m_genMCVecs.m_dSumWeights;
 }
 
+void
+NormIntInterface::applyPoissonWeights( unsigned int seed ) const {
+  if( !m_accMCVecs.m_dataLoaded ) loadMC();
+
+  TRandom3 rndGen(seed);
+
+  vector<int> accW( m_accMCVecs.m_iNTrueEvents );
+  for( size_t i = 0; i < accW.size(); ++i ) accW[i] = rndGen.Poisson(1.0);  
+  m_accMCVecs.applyPoissonWeights( &accW[0] );
+  m_sumAccWeights = m_accMCVecs.m_dSumWeights;
+
+  m_emptyNormIntCache = true;
+}
+
 void
 NormIntInterface::invalidateTerms(){
   
diff --git a/AmpTools/IUAmpTools/NormIntInterface.h b/AmpTools/IUAmpTools/NormIntInterface.h
index 67f38dfc..97872d86 100644
--- a/AmpTools/IUAmpTools/NormIntInterface.h
+++ b/AmpTools/IUAmpTools/NormIntInterface.h
@@ -88,6 +88,7 @@ class NormIntInterface
   // needs to be virtual so parallel implementations can properly
   // override this function
   virtual void forceCacheUpdate( bool normIntOnly = false ) const;
+  virtual void applyPoissonWeights( unsigned int seed ) const;
   
   void invalidateTerms();
 
@@ -103,8 +104,8 @@ class NormIntInterface
   const double* ampIntMatrix() const  { return m_ampIntCache;  }
   const double* normIntMatrix() const { return m_normIntCache; }
   
-  void setGenEvents( double sumWeights ) { m_sumGenWeights = sumWeights; }
-  void setAccEvents( double sumWeights ) { m_sumAccWeights = sumWeights; }
+  void setGenEvents( double sumWeights ) const { m_sumGenWeights = sumWeights; }
+  void setAccEvents( double sumWeights ) const { m_sumAccWeights = sumWeights; }
   
 protected:
   
diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc
index cfcacf8d..a1b27b4f 100644
--- a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc
+++ b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.cc
@@ -276,6 +276,35 @@ LikelihoodCalculatorMPI::computeLikelihood()
   MPI_Send( (double*)&data, 5, MPI_DOUBLE, 0, MPITag::kDoubleSend, MPI_COMM_WORLD );
 }
 
+void
+LikelihoodCalculatorMPI::applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed ){
+  assert( m_isLeader );
+
+  int cmnd[2] = { m_thisId, LikelihoodManagerMPI::kApplyPoissonWeights };
+  MPI_Bcast( cmnd, 2, MPI_INT, 0, MPI_COMM_WORLD );
+
+  unsigned int seeds[2] = { dataSeed, mcSeed };
+  MPI_Bcast( seeds, 2, MPI_UNSIGNED, 0, MPI_COMM_WORLD );
+
+  // despite the leader holding no sig/bkgd data for reweighting, we need to join 
+  // the normIntInterface MPI handshakes to avoid a desync.
+  m_normInt.applyPoissonWeights( mcSeed );
+  
+}
+
+void
+LikelihoodCalculatorMPI::applyBootstrapWeights(){
+  assert( !m_isLeader );
+
+  unsigned int seeds[2];
+  MPI_Bcast( seeds, 2, MPI_UNSIGNED, 0, MPI_COMM_WORLD );
+
+  // give each follower an independent generator
+  unsigned int dataSeed = seeds[0] + 1000003u * (unsigned int) m_rank;
+  unsigned int mcSeed = seeds[1] + 1000003u * (unsigned int) m_rank;
+  LikelihoodCalculator::applyPoissonWeights( dataSeed, mcSeed );
+}
+
 void
 LikelihoodCalculatorMPI::setupMPI()
 {
diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h
index eb88cd6b..c171e920 100644
--- a/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h
+++ b/AmpTools/IUAmpToolsMPI/LikelihoodCalculatorMPI.h
@@ -40,6 +40,7 @@
 #include "IUAmpToolsMPI/MPITag.h"
 #include "IUAmpToolsMPI/ParameterManagerMPI.h"
 #include "IUAmpTools/LikelihoodCalculator.h"
+#include "IUAmpTools/NormIntInterface.h"
 
 class LikelihoodManagerMPI;
 class MISubject;
@@ -137,6 +138,9 @@ class LikelihoodCalculatorMPI : public LikelihoodCalculator
    *  for the entire job summed over all following processes.
    */
   double numSignalEvents();
+
+  virtual void applyPoissonWeights( unsigned int dataSeed, unsigned int mcSeed ); // leader entry point
+  void applyBootstrapWeights(); // follower handler
   
 private:
   
diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc
index 3582a22e..09be58b2 100644
--- a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc
+++ b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.cc
@@ -106,6 +106,11 @@ LikelihoodManagerMPI::deliverLikelihood()
         
         likCalc->computeLikelihood();
         break;
+      
+      case kApplyPoissonWeights:
+
+        likCalc->applyBootstrapWeights();
+        break;
         
       default:
         
diff --git a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h
index 5d1c79ed..4d179739 100644
--- a/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h
+++ b/AmpTools/IUAmpToolsMPI/LikelihoodManagerMPI.h
@@ -60,6 +60,7 @@ class LikelihoodManagerMPI
                     kComputeIntegrals,
                     kUpdateParameters,
                     kUpdateAmpParameter,
+                    kApplyPoissonWeights,
                     kFinalizeFit,
                     kExit };
 
diff --git a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
index 0a5e675e..d230d826 100644
--- a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
+++ b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.cc
@@ -91,17 +91,29 @@ NormIntInterfaceMPI::setupMPI()
   MPI_Comm_rank( MPI_COMM_WORLD, &m_rank );
   MPI_Comm_size( MPI_COMM_WORLD, &m_numProc );
   
-  MPI_Status status;
-  
   m_isLeader = ( m_rank == 0 );
-  
-  double totalGenWeights = 0;
-  double totalAccWeights = 0;
-    
+
+  if( !m_isLeader ){
+    // load the MC into CPU memory on the follower nodes -- this is effectively
+    // a copy from the DataReaderMPI cache to the AmpVecs structure
+    // and it is only run on the followers to avoid large memory
+    // usage on the lead node
+    loadMC();
+  }
+
+  syncEventCounts();
+}
+
+void
+NormIntInterfaceMPI::syncEventCounts() const {
+  MPI_Status status;
+
   if( m_isLeader ){
-    
+
+    double totalGenWeights = 0;
+    double totalAccWeights = 0;
+
     for( int i = 1; i < m_numProc; ++i ){
-      
       double thisWeights;
 
       // trigger sending of events from followers -- data is irrelevant
@@ -124,15 +136,8 @@ NormIntInterfaceMPI::setupMPI()
     
     setGenEvents( totalGenWeights );
     setAccEvents( totalAccWeights );
-  }
+    }
   else{
-    
-    // load the MC into CPU memory on the follower nodes -- this is effectively
-    // a copy from the DataReaderMPI cache to the AmpVecs structure
-    // and it is only run on the followers to avoid large memory
-    // usage on the lead node
-    loadMC();
-
     double thisWeights;
 
     // if we are not the leader, send generated and accepted events
@@ -210,3 +215,9 @@ NormIntInterfaceMPI::loadMC() const{
   if( !m_isLeader ) NormIntInterface::loadMC();
 }
 
+void
+NormIntInterfaceMPI::applyPoissonWeights( unsigned int seed ) const {
+  if( !m_isLeader ) NormIntInterface::applyPoissonWeights( seed );
+  syncEventCounts();
+  forceCacheUpdate( true );
+} 
\ No newline at end of file
diff --git a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.h b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.h
index 067c047f..bec3f269 100644
--- a/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.h
+++ b/AmpTools/IUAmpToolsMPI/NormIntInterfaceMPI.h
@@ -59,11 +59,13 @@ class NormIntInterfaceMPI : public NormIntInterface
   void forceCacheUpdate( bool normIntOnly = false ) const;
 
   void loadMC() const;
+  void applyPoissonWeights( unsigned int seed ) const override;
 
 private:
   
   void setupMPI();
   void sumIntegrals( IntType type ) const;
+  void syncEventCounts() const;
   
   bool m_mpiSetup;
   

From 5ae624bc63d922f52d01d75d4c965187d3a7a7d4 Mon Sep 17 00:00:00 2001
From: kevScheuer <kscheuer@jlab.org>
Date: Fri, 2 Oct 2026 13:16:56 -0400
Subject: [PATCH 2/2] Prevent empty GPU userVars being read

When a persistent AmpToolsInterface was used (e.g.
randomized fitting or bootstrap sampling), the
very first fit would have the correct ampInt, but
all subsequent ampInt's would be bugged.

After the first fit, the userVars on the GPU
appeared to be cleared, and so this commit
re-copies them.
---
 AmpTools/IUAmpTools/AmpVecs.cc | 1 +
 1 file changed, 1 insertion(+)

diff --git a/AmpTools/IUAmpTools/AmpVecs.cc b/AmpTools/IUAmpTools/AmpVecs.cc
index da218428..652fb749 100644
--- a/AmpTools/IUAmpTools/AmpVecs.cc
+++ b/AmpTools/IUAmpTools/AmpVecs.cc
@@ -424,6 +424,7 @@ AmpVecs::allocateTerms( const IntensityManager& intenMan, bool bAllocIntensity,
 #else
   
   m_gpuMan.initTerms( *this, chunkSize );
+  if( m_pdUserVars != 0 && m_userVarsValid ) m_gpuMan.copyUserVarsToGPU( *this );
 
 #endif // GPU_ACCELERATION