diff --git a/src/libraries/AMPTOOLS_AMPS/BreitWigner.cc b/src/libraries/AMPTOOLS_AMPS/BreitWigner.cc
index 3dfc84ef1..aeb0ea75f 100644
--- a/src/libraries/AMPTOOLS_AMPS/BreitWigner.cc
+++ b/src/libraries/AMPTOOLS_AMPS/BreitWigner.cc
@@ -1,5 +1,3 @@
-
-
 #include <cassert>
 #include <iostream>
 #include <string>
diff --git a/src/libraries/AMPTOOLS_AMPS/GPUIso_ps_refl_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPUIso_ps_refl_kernel.cu
new file mode 100644
index 000000000..f4578a11f
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/GPUIso_ps_refl_kernel.cu
@@ -0,0 +1,18 @@
+#include <stdio.h>
+
+#include "GPUManager/GPUCustomTypes.h"
+#include "GPUManager/CUDA-Complex.cuh"
+
+__global__ void GPUIso_ps_refl_kernel( GPU_AMP_PROTO )
+{
+	int iEvent = GPU_THIS_EVENT;
+
+	WCUComplex ans = { GPU_UVARS(0), GPU_UVARS(1) };
+    	pcDevAmp[iEvent] = ans;	
+}
+
+void GPUIso_ps_refl_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO )
+{  
+  	GPUIso_ps_refl_kernel<<< dimGrid, dimBlock >>>( GPU_AMP_ARGS );
+}
+
diff --git a/src/libraries/AMPTOOLS_AMPS/GPUPiPiSWaveAMPK_kernel.cu b/src/libraries/AMPTOOLS_AMPS/GPUPiPiSWaveAMPK_kernel.cu
new file mode 100644
index 000000000..ab8f1bb8e
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/GPUPiPiSWaveAMPK_kernel.cu
@@ -0,0 +1,19 @@
+#include <stdio.h>
+
+#include "GPUManager/GPUCustomTypes.h"
+#include "GPUManager/CUDA-Complex.cuh"
+
+__global__ void GPUPiPiSWaveAMPK_kernel( GPU_AMP_PROTO )
+{
+	int iEvent = GPU_THIS_EVENT;
+
+	WCUComplex ans = { GPU_UVARS(0), GPU_UVARS(1) };
+    	pcDevAmp[iEvent] = ans;	
+}
+
+void GPUPiPiSWaveAMPK_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO )
+{  
+
+  	GPUPiPiSWaveAMPK_kernel<<< dimGrid, dimBlock >>>( GPU_AMP_ARGS );
+
+}
\ No newline at end of file
diff --git a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc
new file mode 100644
index 000000000..f83c42af7
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.cc
@@ -0,0 +1,234 @@
+
+#include <cassert>
+#include <iostream>
+#include <string>
+#include <sstream>
+#include <cstdlib>
+
+#include "TLorentzVector.h"
+#include "TLorentzRotation.h"
+#include "TFile.h"
+
+#include "IUAmpTools/Kinematics.h"
+#include "AMPTOOLS_AMPS/Iso_ps_refl.h"
+#include "AMPTOOLS_AMPS/clebschGordan.h"
+#include "AMPTOOLS_AMPS/wignerD.h"
+#include "AMPTOOLS_AMPS/decayAngles.h"
+#include "AMPTOOLS_AMPS/barrierFactor.h"
+
+#include "UTILITIES/BeamProperties.h"
+
+// Function to check if a string is a valid number
+static bool isValidNumber(const string& argInput, double &value){
+    char* end = nullptr;
+    errno = 0;  // reset global error
+    value = std::strtod(argInput.c_str(), &end);
+
+    // Check if 
+    // (1) no conversion was performed 
+    // (2) there are leftover characters
+    // (3) an overflow/underflow occurred   
+    if(end == argInput.c_str() || *end != '\0' || errno != 0) {
+        return false;  // not a valid number
+    }
+    // If end points to the end of string, it's fully numeric
+    return true;
+}
+
+static double parseValidatedNumber(const string& label, const string& argInput){
+    double tmpValue = 0.0;
+    if(!isValidNumber(argInput, tmpValue)){
+      throw std::invalid_argument("Iso_ps_refl: invalid " + label + ": " + argInput);
+    }
+    return tmpValue;
+}
+
+
+Iso_ps_refl::Iso_ps_refl( const vector< string >& args ) :
+UserAmplitude< Iso_ps_refl >( args ){
+
+  // This function is only for the two-body Isobar decay: xi -> pipi
+  
+  // 6 possibilities to initialize this amplitude:
+  // <S>: isobar spin
+  // <J>: total spin, 
+  // <m>: total spin projection, 
+  // <l>: partial wave, 
+  // <r>: +1/-1 for real/imaginary part; 
+  // <s>: +1/-1 sign in P_gamma term
+
+  // Isobar spin S
+  m_s = static_cast<int>(parseValidatedNumber("S",args[0]));    
+  // Resonance spin J
+  m_j = static_cast<int>(parseValidatedNumber("J", args[1])); 
+  // Resonance spin projection 
+  m_m = static_cast<int>(parseValidatedNumber("M", args[2])); 
+  // Partial wave L
+  m_l = static_cast<int>(parseValidatedNumber("L", args[3])); 
+  // Real (+1) or imaginary (-1)
+  m_real = static_cast<int>(parseValidatedNumber("Re/Im", args[4])); 
+  // Sign for polarization in amplitude
+  m_sign = static_cast<int>(parseValidatedNumber("P_gamma sign", args[5])); 
+
+  
+  // make sure values are reasonable
+  assert(m_s >= 0 && m_s <= 2);
+  assert( abs( m_m ) <= m_j );
+  // m_real = +1 for real
+  // m_real = -1 for imag
+  assert( abs( m_real ) == 1 );
+  // m_sign = +1 for 1 + Pgamma
+  // m_sign = -1 for 1 - Pgamma
+  assert( abs( m_sign ) == 1 );
+
+  // Default polarization information stored in tree
+  m_polInTree = true;
+
+  // Loop over any additional amplitude arguments to change defaults
+  for(uint ioption=6; ioption<args.size(); ioption++) {
+	  TString option = args[ioption].c_str();
+    // Polarization provided in configuration file
+    if(ioption==6){
+      m_polInTree = false;
+
+      polAngle = parseValidatedNumber("polarization angle", args[6]);    
+
+      TString polOption = args[7].c_str();
+      if(polOption.Contains(".root")){
+        polFraction = 0.0;
+        TFile* f = new TFile(polOption);
+        polFrac_vs_E = (TH1D*)f->Get(args[8].c_str());
+        if(polFrac_vs_E  != nullptr ){
+          throw std::runtime_error(
+            "Iso_ps_refl ERROR: Could not find histogram '" + args[8] +
+            "' in file " + std::string(polOption.Data()));
+        }
+      }
+      else{
+        polFraction = parseValidatedNumber("polarization fraction", args[7]);
+      }
+    }
+  }  
+}
+
+void Iso_ps_refl::calcUserVars( GDouble** pKin, GDouble* userVars ) const{
+
+  TLorentzVector beam;
+  TVector3 eps;
+  GDouble beam_polFraction;
+  GDouble beam_polAngle;
+
+  if(m_polInTree){
+    beam.SetPxPyPzE( 0.0, 0.0, pKin[0][0], pKin[0][0]);
+    eps.SetXYZ(pKin[0][1], pKin[0][2], 0.0); // beam polarization vector;
+
+    beam_polFraction = eps.Mag();
+    beam_polAngle = eps.Phi();
+  }
+  else{
+    beam.SetPxPyPzE( pKin[0][1], pKin[0][2], pKin[0][3], pKin[0][0] );
+    beam_polAngle = polAngle;
+
+    if(polFraction > 0.0){ // for fitting with fixed polarization
+	    beam_polFraction = polFraction;
+    }
+    else{ // for fitting with polarization vs E_gamma from input histogram 
+	    int bin = polFrac_vs_E->GetXaxis()->FindBin(pKin[0][0]);
+	    if (bin == 0 || bin > polFrac_vs_E->GetXaxis()->GetNbins()){
+		    beam_polFraction = 0.0;
+	    } else
+		    beam_polFraction = polFrac_vs_E->GetBinContent(bin);
+    }
+  }
+
+
+  // Fill in four-vectors for final state particles
+  TLorentzVector target(0,0,0,0.938272);   // Fixed target
+  TLorentzVector recoil( pKin[1][1]+pKin[5][1], pKin[1][2]+pKin[5][2], pKin[1][3]+pKin[5][3], pKin[1][0]+pKin[5][0] ); // Recoiling baryon
+  TLorentzVector ps(pKin[2][1], pKin[2][2], pKin[2][3], pKin[2][0]);     // 1st after proton is always the pseudoscalar meson
+
+  
+  // Compute isobar meson from its decay products (xi -> pipi)
+  // Make sure the order of daughters is correct in the config file!
+  TLorentzVector iso_daught1(pKin[3][1], pKin[3][2], pKin[3][3], pKin[3][0]);
+  TLorentzVector iso_daught2(pKin[4][1], pKin[4][2], pKin[4][3], pKin[4][0]);
+  TLorentzVector iso = iso_daught1 + iso_daught2;
+
+  // Final meson system P4
+  TLorentzVector X = iso + ps;
+
+
+  //Calculate production angle in the Gottfried-Jackson frame
+  GDouble prod_angle = getPhiProd(beam_polAngle, X, beam, target, 2, true);
+
+  // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame  
+  vector <double> thetaPhiAnglesTwoStep;
+  thetaPhiAnglesTwoStep = getTwoStepAngles(X, iso, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true);
+
+  
+  
+  GDouble cosTheta = TMath::Cos(thetaPhiAnglesTwoStep[0]);
+  GDouble phi = thetaPhiAnglesTwoStep[1];
+  GDouble cosThetaH = TMath::Cos(thetaPhiAnglesTwoStep[2]);
+  GDouble phiH = thetaPhiAnglesTwoStep[3];
+  GDouble m_X = X.M();
+  GDouble m_iso = iso.M();
+  GDouble m_ps = ps.M();
+
+  complex <GDouble> amplitude(0,0);
+  complex <GDouble> i(0,1);
+
+  for (int lambda = -m_s; lambda <= m_s; lambda++) { // sum over helicities
+	  GDouble hel_amp = clebschGordan(m_l, m_s, 0, lambda, m_j, lambda);
+          amplitude += conj(wignerD(m_j, m_m, lambda, cosTheta, phi))*hel_amp*conj(wignerD(m_s, lambda, 0, cosThetaH, phiH));
+  }
+
+  
+  GDouble factor = sqrt(1 + m_sign * beam_polFraction);
+  complex <GDouble> zjm = 0;
+  // - -> + in prod_angle
+  complex <GDouble> rotateY = polar((GDouble)1., (GDouble)(-1. * prod_angle ));  
+
+  if (m_real == 1)
+	  zjm = real(amplitude * rotateY);
+  if (m_real == -1) 
+	  zjm = i*imag(amplitude * rotateY);
+
+  // E852 Nozar thesis has sqrt(2*s+1)*sqrt(2*l+1)*F_l(p_omega)*sqrt(omega)
+  double kinFactor = barrierFactor(m_X, m_l, m_iso, m_ps);
+  //kinFactor *= sqrt(2.*m_s + 1.) * sqrt(2.*m_l + 1.);
+  factor *= kinFactor;
+
+  userVars[uv_ampRe] = ( factor * zjm ).real();
+  userVars[uv_ampIm] = ( factor * zjm ).imag();
+
+  return;
+}
+
+
+/////////////////////// Amplitude Calculation //////////////////////////
+
+complex< GDouble >
+Iso_ps_refl::calcAmplitude( GDouble** pKin, GDouble* userVars ) const
+{
+  return complex< GDouble >( userVars[uv_ampRe], userVars[uv_ampIm] );
+}
+
+void Iso_ps_refl::updatePar( const AmpParameter& par ){
+
+  // could do expensive calculations here on parameter updates  
+}
+
+
+#ifdef GPU_ACCELERATION
+
+void
+Iso_ps_refl::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const {
+
+	GPUIso_ps_refl_exec( dimGrid, dimBlock, GPU_AMP_ARGS );
+
+}
+
+#endif
+
+
diff --git a/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.h b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.h
new file mode 100644
index 000000000..e67be10a6
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/Iso_ps_refl.h
@@ -0,0 +1,93 @@
+#if !defined(ISO_PS_REFL)
+#define ISO_PS_REFL
+
+#include "IUAmpTools/Amplitude.h"
+#include "IUAmpTools/UserAmplitude.h"
+#include "IUAmpTools/AmpParameter.h"
+#include "GPUManager/GPUCustomTypes.h"
+
+#include "TH1D.h"
+#include <string>
+#include <complex>
+#include <vector>
+
+using std::complex;
+using namespace std;
+
+#ifdef GPU_ACCELERATION
+void
+GPUIso_ps_refl_exec( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO );
+#endif
+
+class Kinematics;
+
+class Iso_ps_refl : public UserAmplitude< Iso_ps_refl >
+{
+    
+public:
+	
+	Iso_ps_refl() : UserAmplitude< Iso_ps_refl >() { };
+	Iso_ps_refl( const vector< string >& args );
+        Iso_ps_refl( int m_s, int m_j, int m_m, int m_l, int m_real, int m_sign );
+	
+	string name() const { return "Iso_ps_refl"; }
+    
+	complex< GDouble > calcAmplitude( GDouble** pKin, GDouble* userVars ) const;
+
+	// **********************
+	// The following lines are optional and can be used to precalculate
+	// user-defined data that the amplitudes depend on.
+	
+	// Use this for indexing a user-defined data array and notifying
+	// the framework of the number of user-defined variables.
+		
+	//enum UserVars { uv_cosTheta = 0, uv_Phi, uv_cosThetaH, uv_PhiH,
+	//                uv_prod_Phi, uv_MX, uv_MVec, uv_MPs, uv_beam_polFraction,
+	//                uv_beam_polAngle, kNumUserVars };
+	enum UserVars { uv_ampRe = 0, uv_ampIm, kNumUserVars };
+	unsigned int numUserVars() const { return kNumUserVars; }
+	
+	// This function needs to be defined -- see comments and discussion
+	// in the .cc file.
+	void calcUserVars( GDouble** pKin, GDouble* userVars ) const;
+	
+	// This is an optional addition if the calcAmplitude routine
+	// can run with only the user-defined data and not the original
+	// four-vectors.  It is used to optimize memory usage in GPU
+	// based fits.
+	bool needsUserVarsOnly() const { return true; }
+	
+	// This is an optional addition if the UserVars are the same for each 	
+	// instance of an amplitude.  If it is not used, the memory footprint
+	// grows dramatically as UserVars values are stored for each instance
+	// of the amplitude.  NOTE: To use this make sure that UserVars only 
+	// depend on kinematics and no arguments provided to the amplitude!
+	bool areUserVarsStatic() const { return false; }
+
+	void updatePar( const AmpParameter& par );
+
+#ifdef GPU_ACCELERATION
+
+	void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const;
+
+        bool isGPUEnabled() const { return true; }
+
+#endif // GPU_ACCELERATION
+	
+private:
+        
+        int m_s;
+        int m_j;
+	int m_m;
+	int m_l;
+	int m_real;
+	int m_sign;
+  
+	//AmpParameter polAngle;
+	GDouble polFraction;
+	GDouble polAngle;
+	bool m_polInTree;
+	TH1D *polFrac_vs_E;
+};
+
+#endif
diff --git a/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.cc b/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.cc
new file mode 100644
index 000000000..947195b23
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.cc
@@ -0,0 +1,132 @@
+#include <cassert>
+#include <cstdlib>
+#include <iostream>
+#include <string>
+#include <complex>
+
+#include "TLorentzVector.h"
+#include "IUAmpTools/Kinematics.h"
+#include "breakupMomentumComplex.h"
+#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h"
+
+
+// S-wave dipion mass parametrization amplitude based on:
+//   K.L. Au, D. Morgan and M.R. Pennington," Meson dynamics beyond the quark model: Study of final-state interactions",
+//   PRD 35 1987
+//   with the modification by I. Kachaev from [arxiv:2305.11711] 
+//   Usage:
+//   amplitude <reaction>::<sum>::<ampName> PiPiSWaveAMPK <id_daughter1> <id_daughter2> 
+
+using namespace std;
+
+PiPiSWaveAMPK::PiPiSWaveAMPK( const vector <string> &args ) : UserAmplitude<PiPiSWaveAMPK>(args)
+{
+
+   assert( args.size() == 2 );
+   m_daughters = pair <string,string> (args[0],args[1]);
+
+   setParametrization();
+
+   // need to register a free parameters so the framework knows about it
+   // registerParameter(m_mass);
+
+}
+
+
+void PiPiSWaveAMPK::setParametrization() 
+{
+  
+  //Adler pole and residue in it
+  s0 = -0.0074;  // AMP Table 1, M solution: s_0
+  a11 = 0.1131; // AMP Table 1, M solution: f_2^2
+  
+  //polynom up to the 3rd order
+  c11.resize(4);
+
+  c11[0] = 0.0337;                 // AMP Table 1, M solution: c_11^0
+  c11[1] = -0.3185;                // AMP Table 1, M solution: c_11^1
+  c11[2] = -0.0942;                // AMP Table 1, M solution: c_11^2
+  c11[3] = -0.5927;                // AMP Table 1, M solution: c_11^3
+
+}
+
+
+
+/////////////////////// Amplitude Calculation //////////////////////////
+
+complex <GDouble> PiPiSWaveAMPK::calcAmplitude( GDouble** pKin ) const
+{
+  TLorentzVector pion1_P4, pion2_P4, dipion_P4, temp_P4;
+  
+  for( unsigned int ii = 0; ii < m_daughters.first.size(); ++ii ){
+
+    char num[2]= {m_daughters.first[ii], '\0'};
+    int index = atoi(num);
+
+    temp_P4.SetPxPyPzE(pKin[index][1], pKin[index][2], pKin[index][3], pKin[index][0]);
+    pion1_P4 += temp_P4;
+    dipion_P4 += temp_P4;
+  }
+  
+  for( unsigned int ii = 0; ii < m_daughters.second.size(); ++ii ){
+    
+    char num[2]= {m_daughters.second[ii], '\0'};
+    int index = atoi(num);
+
+    temp_P4.SetPxPyPzE(pKin[index][1], pKin[index][2], pKin[index][3], pKin[index][0]);
+    pion2_P4 += temp_P4;
+    dipion_P4 += temp_P4;
+  }
+
+
+  GDouble mass  = dipion_P4.M();
+  GDouble s  = dipion_P4.M2();
+   
+
+  if (fabs(s - s0) < 1.e-6){
+     mass += 1.e-6;
+     s = mass*mass;
+  }
+
+  
+  const complex <GDouble> qPiPi   = breakupMomentumComplex(mass, _piChargedMass,   _piChargedMass );
+  const complex <GDouble> qPi0Pi0 = breakupMomentumComplex(mass, _piNeutralMass,   _piNeutralMass );
+  const GDouble scale = (s / (4*_kaonMeanMass*_kaonMeanMass)) - 1;
+
+
+  complex <GDouble> rho11 = (qPiPi + qPi0Pi0) / mass;
+
+  complex <GDouble> imag(0.,1.);  
+  complex <GDouble> M11(0.,0.),T11(0.,0.);
+
+  //Add the Adler pole term
+   M11 += a11/(s - s0);
+
+  //Add the polynomial terms  
+  for (unsigned int ii = 0; ii < c11.size(); ++ii) 
+    M11 += c11[ii]*pow(scale, (int)ii);
+
+  //Now calculate the T_{11} matrix element using the T11 to M11 relation
+  T11 = 1./(M11 - imag*rho11); 
+
+
+  return T11;
+}
+
+
+
+
+//This function may be used instead of the separate amplitude ' breakupMomentumComplex' 
+template<typename mType> complex<GDouble> PiPiSWaveAMPK::breakupMom(mType m, GDouble mDec1, GDouble mDec2) const{
+   complex<GDouble> result = PiPiSWaveAMPK::phaseSpaceFac(m, mDec1, mDec2)*m/GDouble(2.);
+   return result;  
+}
+
+
+#ifdef GPU_ACCELERATION
+void PiPiSWaveAMPK::launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const {
+
+  GPUPiPiSWaveAMPK_exec( dimGrid,  dimBlock, GPU_AMP_ARGS);
+  
+}
+#endif 
diff --git a/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.h b/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.h
new file mode 100644
index 000000000..3a813e276
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/PiPiSWaveAMPK.h
@@ -0,0 +1,71 @@
+#if !defined(PIPISWAVEAMPK)
+#define PIPISWAVEAMPK
+
+#include "particleType.h"
+#include "IUAmpTools/UserAmplitude.h"
+#include "IUAmpTools/AmpParameter.h"
+#include "IUAmpTools/Amplitude.h"
+#include "GPUManager/GPUCustomTypes.h"
+
+#include <utility>
+#include <string>
+#include <complex>
+#include <vector>
+#include <iostream>
+
+
+using std::complex;
+using namespace std;
+
+
+#ifdef GPU_ACCELERATION
+void   GPUPiPiSWaveAMPK_exec(dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO );
+#endif 
+
+
+
+class Kinematics;
+
+class PiPiSWaveAMPK : public UserAmplitude<PiPiSWaveAMPK>{
+
+   public:
+
+      PiPiSWaveAMPK() : UserAmplitude <PiPiSWaveAMPK>() { }
+      PiPiSWaveAMPK( const vector<string> &args );
+      ~PiPiSWaveAMPK(){}
+
+      string name() const { return "PiPiSWaveAMPK"; }
+      void setParametrization(); 
+      complex <GDouble> calcAmplitude( GDouble** pKin ) const;
+
+#ifdef GPU_ACCELERATION
+
+     void launchGPUKernel( dim3 dimGrid, dim3 dimBlock, GPU_AMP_PROTO ) const;
+     bool isGPUEnabled() const { return true; }
+
+#endif 
+
+   private:
+
+      pair< string, string > m_daughters;    
+
+      GDouble s0; 
+      GDouble a11; 
+      std::vector <GDouble> c11; 
+
+  
+      const GDouble _piChargedMass = ParticleMass(PiPlus); //0.13957039;
+      const GDouble _piNeutralMass = ParticleMass(Pi0); //0.13497680;
+      const GDouble _kaonChargedMass = ParticleMass(KPlus); //0.49367700;
+      const GDouble _kaonNeutralMass = 0.5*(ParticleMass(KLong) + ParticleMass(KShort)); //0.49761400;
+      const GDouble _kaonMeanMass = 0.5*(_kaonChargedMass + _kaonNeutralMass); //0.5*(0.49367700 + 0.49761400);
+
+
+  
+      complex<GDouble> phaseSpaceFac( GDouble m, GDouble mDec1, GDouble mDec2 ) const;
+      template<typename mType>
+      complex<GDouble> breakupMom( mType m, GDouble mDec1, GDouble mDec2 ) const;
+
+};
+
+#endif
diff --git a/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.cc b/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.cc
new file mode 100644
index 000000000..afe441013
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.cc
@@ -0,0 +1,20 @@
+#include <cmath>
+#include <complex>
+#include "breakupMomentumComplex.h"
+
+// mass0 = mass of parent
+// mass1 = mass of first daughter
+// mass2 = mass of second daughter
+
+complex <GDouble> breakupMomentumComplex( GDouble mass0, GDouble mass1, GDouble mass2 ){
+	
+  complex <GDouble> q = std::sqrt(    mass0*mass0*mass0*mass0 + 
+						  mass1*mass1*mass1*mass1 +
+						  mass2*mass2*mass2*mass2 -
+						  2.0*mass0*mass0*mass1*mass1 -
+						  2.0*mass0*mass0*mass2*mass2 -
+						  2.0*mass1*mass1*mass2*mass2  ) / (2.0 * mass0);
+
+   return q;
+	
+}
diff --git a/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.h b/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.h
new file mode 100644
index 000000000..c80bcb451
--- /dev/null
+++ b/src/libraries/AMPTOOLS_AMPS/breakupMomentumComplex.h
@@ -0,0 +1,17 @@
+#include "IUAmpTools/Amplitude.h"
+#include "IUAmpTools/AmpParameter.h"
+#include "IUAmpTools/UserAmplitude.h"
+#include "GPUManager/GPUCustomTypes.h"
+
+
+
+#if !defined(BREAKUPMOMENTUMCOMPLEX)
+#define BREAKUPMOMENTUMCOMPLEX
+
+// mass0 = mass of parent
+// mass1 = mass of first daughter
+// mass2 = mass of second daughter
+
+complex <GDouble> breakupMomentumComplex( GDouble mass0, GDouble mass1, GDouble mass2 );
+
+#endif
diff --git a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc
new file mode 100644
index 000000000..956c633fc
--- /dev/null
+++ b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.cc
@@ -0,0 +1,207 @@
+#include "TLorentzVector.h"
+#include "TLorentzRotation.h"
+
+#include "AMPTOOLS_AMPS/decayAngles.h"
+
+#include "AMPTOOLS_DATAIO/IsoPsPlotGenerator.h"
+#include "IUAmpTools/Histogram1D.h"
+#include "IUAmpTools/Kinematics.h"
+
+
+
+// Function to check if a string is a valid number
+static bool isValidNumber(const string& argInput, double &value){
+    char* end = nullptr;
+    errno = 0;  // reset global error
+    value = std::strtod(argInput.c_str(), &end);
+
+    // Check if 
+    // (1) no conversion was performed 
+    // (2) there are leftover characters
+    // (3) an overflow/underflow occurred   
+    if(end == argInput.c_str() || *end != '\0' || errno != 0) {
+        return false;  // not a valid number
+    }
+    // If end points to the end of string, it's fully numeric
+    return true;
+}
+
+
+static double parseValidatedNumber(const string& label, const string& argInput){
+    double tmpValue = 0.0;
+    if(!isValidNumber(argInput, tmpValue)){
+      throw std::invalid_argument("Iso_ps_refl: invalid " + label + ": " + argInput);
+    }
+    return tmpValue;
+}
+
+
+
+
+/* Constructor to display FitResults */
+IsoPsPlotGenerator::IsoPsPlotGenerator( const FitResults& results, Option opt ) :
+PlotGenerator( results, opt )
+{
+	createHistograms();
+}
+
+IsoPsPlotGenerator::IsoPsPlotGenerator( const FitResults& results ) :
+PlotGenerator( results )
+{
+	createHistograms();
+}
+
+
+/* Constructor for event generator (no FitResult) */
+IsoPsPlotGenerator::IsoPsPlotGenerator( ) :
+PlotGenerator( )
+{
+	createHistograms();
+}
+
+void IsoPsPlotGenerator::createHistograms( ) {
+  cout << " calls to bookHistogram go here" << endl;
+  
+   bookHistogram( kProd_Ang, new Histogram1D( 50, -180., 180., "ProdAng", "Production Angle [deg.]" ) );
+   bookHistogram( kCosTheta, new Histogram1D( 50, -1., 1., "CosTheta_GJ", "cos#theta^{[GJ]}" ) );
+   bookHistogram( kPhi, new Histogram1D( 50, -180., 180., "Phi_GJ", "#phi^{[GJ]} [deg.]" ) );
+   bookHistogram( kCosThetaH, new Histogram1D( 50, -1., 1., "CosTheta_HF", "cos#theta^{[HF]}" ) );
+   bookHistogram( kPhiH, new Histogram1D( 50, -180., 180., "Phi_HF", "#phi^{[HF]} [deg.]" ) );
+
+   bookHistogram( kIsoMass, new Histogram1D( 200, 0., 3., "MIso", "m(2#pi)  [GeV]") );
+   bookHistogram( kIsoPsMass, new Histogram1D( 200, 0.2, 3.2, "MIsoPs", "m(3#pi)  [GeV]") );
+   bookHistogram( kt, new Histogram1D( 200, 0, 1.0 , "t", "-t  [GeV^{2}]" ) );
+   bookHistogram( kRecoilMass, new Histogram1D( 200, 0.9, 1.9 , "ProtonPiplusL_M", "m(p#pi^{+}_{L}) [GeV]" ) );
+   
+   bookHistogram( kProtonPsMass, new Histogram1D( 200, 0.8, 3.8, "ProtonPiminus_M", "m(p#pi^{-}) [GeV]" ) );
+   bookHistogram( kRecoilPsMass, new Histogram1D( 200, 1.0, 4.0, "ProtonPiplusLPiminus_M", "m(p#pi^{+}_{L}#pi^{-}) [GeV]" ) );
+ 
+}
+
+void
+IsoPsPlotGenerator::projectEvent( Kinematics* kin ){
+
+  // this function will make this class backwards-compatible with older versions
+  // (v0.10.x and prior) of AmpTools, but will not be able to properly obtain
+  // the polariation plane in the lab when multiple orientations are used
+  projectEvent( kin, "" );
+}
+
+void
+IsoPsPlotGenerator::projectEvent( Kinematics* kin, const string& reactionName ){
+
+   // We work only with a 2-body vector decay 
+
+  
+  // Fixed target
+   TLorentzVector target(0,0,0,0.938272);
+
+   
+   //cout << "project event" << endl;
+   TLorentzVector beam   = kin->particle( 0 );
+   TLorentzVector proton = kin->particle( 1 );
+   TLorentzVector bach = kin->particle( 2 );
+   TLorentzVector iso_daught1 = kin->particle( 3 );
+   TLorentzVector iso_daught2 = kin->particle( 4 );
+   TLorentzVector piplusL = kin->particle( 5 );
+   
+
+   // Final state P4 momenta
+   TLorentzVector X = iso_daught1 + iso_daught2 + bach;
+   TLorentzVector recoil = proton + piplusL;
+
+   //Momenta for the 1st permutation   
+   TLorentzVector iso_a = iso_daught1 + iso_daught2;
+   TLorentzVector proton_ps_a = proton + bach;
+   TLorentzVector recoil_ps_a = recoil + bach;
+   
+   //Momenta for the 2nd permutation (bach <=> iso_daught1)  
+   TLorentzVector iso_b = bach + iso_daught2;
+   TLorentzVector proton_ps_b = proton + iso_daught1;
+   TLorentzVector recoil_ps_b = recoil + iso_daught1;
+
+   
+
+   // Properly read polarization angle from config file if provided
+   double beam_polAngle=0;
+   // Check config file for optional parameters -- we assume here that the first amplitude in the list is a Iso_ps_refl amplitude
+   const vector< string > args = cfgInfo()->amplitudeList( reactionName, "", "" ).at(0)->factors().at(0);
+   for(uint ioption=5; ioption<args.size(); ioption++) {
+          TString option = args[ioption].c_str();
+	  if(ioption == 6) beam_polAngle = parseValidatedNumber("polarization angle", args[6]);
+   }
+
+
+
+   //Momentum transfer
+   double Mandt = fabs((target-recoil).M2());
+   
+   //Calculate production angle in the Gottfried-Jackson frame
+   double prod_angle = TMath::RadToDeg()*getPhiProd(beam_polAngle, X, beam, target, 2, true);
+
+   // Calculate decay angles for X in the Gottfried-Jackson frame and for Isobar in the Helicity frame  
+   // Angles for the 1st permutation
+   vector <double> thetaPhiAnglesTwoStep_a = getTwoStepAngles(X, iso_a, iso_daught1, TLorentzVector(0,0,0,0), beam, target, 2, true);
+
+   // Angles for the 2nd permutation (bach <=> iso_daught1)
+   vector <double> thetaPhiAnglesTwoStep_b = getTwoStepAngles(X, iso_b, bach, TLorentzVector(0,0,0,0), beam, target, 2, true);
+
+   
+   //Symmetrized angles and masses will be passed as vectors of unit length   
+   vector <double> cosTheta_a = {TMath::Cos(thetaPhiAnglesTwoStep_a[0])};
+   vector <double> cosTheta_b = {TMath::Cos(thetaPhiAnglesTwoStep_b[0])};
+
+   vector <double> phi_a = {TMath::RadToDeg()*thetaPhiAnglesTwoStep_a[1]};
+   vector <double> phi_b = {TMath::RadToDeg()*thetaPhiAnglesTwoStep_b[1]};
+
+   vector <double> cosThetaH_a = {TMath::Cos(thetaPhiAnglesTwoStep_a[2])};
+   vector <double> cosThetaH_b = {TMath::Cos(thetaPhiAnglesTwoStep_b[2])};
+   
+   vector <double> phiH_a = {TMath::RadToDeg()*thetaPhiAnglesTwoStep_a[3]};
+   vector <double> phiH_b = {TMath::RadToDeg()*thetaPhiAnglesTwoStep_b[3]};
+   
+   vector <double> iso_mass_a = {iso_a.M()};
+   vector <double> iso_mass_b = {iso_b.M()};
+   
+   vector <double> protonps_mass_a = {proton_ps_a.M()};
+   vector <double> protonps_mass_b = {proton_ps_b.M()};
+   
+   vector <double> recoilps_mass_a = {recoil_ps_a.M()};
+   vector <double> recoilps_mass_b = {recoil_ps_b.M()};
+
+   
+
+   //First, insensitive to the symmetrization quantities
+   fillHistogram( kProd_Ang, prod_angle );
+   fillHistogram( kt, Mandt );
+   fillHistogram( kRecoilMass, recoil.M() );
+   fillHistogram( kIsoPsMass, X.M() );
+
+
+   //Now, symmetrized quantities 
+   fillHistogram( kCosTheta, cosTheta_a, 0.5 );
+   fillHistogram( kCosTheta, cosTheta_b, 0.5 );
+ 
+   fillHistogram( kPhi, phi_a, 0.5 );
+   fillHistogram( kPhi, phi_b, 0.5 );
+ 
+   fillHistogram( kCosThetaH, cosThetaH_a, 0.5 );
+   fillHistogram( kCosThetaH, cosThetaH_b, 0.5 );
+ 
+   fillHistogram( kPhiH, phiH_a, 0.5 );
+   fillHistogram( kPhiH, phiH_b, 0.5 );
+
+   
+   fillHistogram( kIsoMass, iso_mass_a, 0.5 );
+   fillHistogram( kIsoMass, iso_mass_b, 0.5 );
+  
+   fillHistogram( kProtonPsMass, protonps_mass_a, 0.5 );
+   fillHistogram( kProtonPsMass, protonps_mass_b, 0.5 );
+  
+   fillHistogram( kRecoilPsMass, recoilps_mass_a, 0.5 );
+   fillHistogram( kRecoilPsMass, recoilps_mass_b, 0.5 );
+
+
+   
+   
+}
diff --git a/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.h b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.h
new file mode 100644
index 000000000..0f5123bad
--- /dev/null
+++ b/src/libraries/AMPTOOLS_DATAIO/IsoPsPlotGenerator.h
@@ -0,0 +1,36 @@
+#if !(defined ISOPSPLOTGENERATOR)
+#define ISOPSPLOTGENERATOR
+
+#include <vector>
+#include <string>
+
+#include "IUAmpTools/PlotGenerator.h"
+
+using namespace std;
+
+class FitResults;
+class Kinematics;
+
+class IsoPsPlotGenerator : public PlotGenerator
+{
+    
+public:
+  
+  // create an index for different histograms
+  enum { kProd_Ang = 0, kCosTheta = 1, kPhi = 2, kCosThetaH = 3, kPhiH = 4, kIsoMass = 5, kIsoPsMass = 6, kt = 7, kRecoilMass = 8, kProtonPsMass = 9, kRecoilPsMass = 10, kNumHists};
+
+  IsoPsPlotGenerator( const FitResults& results, Option opt);
+  IsoPsPlotGenerator( const FitResults& results );
+  IsoPsPlotGenerator( );
+ 
+private:
+  
+  void projectEvent( Kinematics* kin );
+  void projectEvent( Kinematics* kin, const string& reactionName );
+
+  void createHistograms( );
+ 
+};
+
+#endif
+
diff --git a/src/programs/AmplitudeAnalysis/SConscript b/src/programs/AmplitudeAnalysis/SConscript
index 2e64c9056..d2149bdff 100644
--- a/src/programs/AmplitudeAnalysis/SConscript
+++ b/src/programs/AmplitudeAnalysis/SConscript
@@ -3,10 +3,11 @@ import sbms
 
 Import('*')
 
-subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', 'twolepton_plotter', 'twoleptonGJ_plotter', 'split_mass', 'split_t', 'threepi_plotter_schilling', 'omega_radiative_plotter', 'project_moments', 'plot_etapi_delta', 'project_moments_polarized', 'Bootstrap_plot_etapi_delta_SPDG_allamps_mass_t_bins', 'Pol_moments_viafittedPW', 'project_moments_SPD_etapi0_posepsilon', 'omegapi_plotter', 'vecps_plotter', 'plot_etapi0'] 
+subdirs = ['fit', 'twopi_plotter', 'twopi_plotter_amp', 'twopi_plotter_mom', 'twopi_plotter_primakoff', 'twolepton_plotter', 'twoleptonGJ_plotter', 'split_mass', 'split_t', 'threepi_plotter_schilling', 'omega_radiative_plotter', 'project_moments', 'plot_etapi_delta', 'project_moments_polarized', 'Bootstrap_plot_etapi_delta_SPDG_allamps_mass_t_bins', 'Pol_moments_viafittedPW', 'project_moments_SPD_etapi0_posepsilon', 'omegapi_plotter', 'vecps_plotter', 'vecps3pi_plotter', 'isops_plotter', 'plot_etapi0'] 
 
 SConscript(dirs=subdirs, exports='env osname', duplicate=0)
 
+#fitMPI
 # Optional targets
 optdirs = ['fitMPI']
 sbms.OptionallyBuild(env, optdirs)
diff --git a/src/programs/AmplitudeAnalysis/fit/fit.cc b/src/programs/AmplitudeAnalysis/fit/fit.cc
index c630f946d..b723de23d 100644
--- a/src/programs/AmplitudeAnalysis/fit/fit.cc
+++ b/src/programs/AmplitudeAnalysis/fit/fit.cc
@@ -37,6 +37,7 @@
 #include "AMPTOOLS_AMPS/Zlm.h"
 #include "AMPTOOLS_AMPS/BreitWigner.h"
 #include "AMPTOOLS_AMPS/BreitWigner3body.h"
+#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h"
 #include "AMPTOOLS_AMPS/b1piAngAmp.h"
 #include "AMPTOOLS_AMPS/Uniform.h"
 #include "AMPTOOLS_AMPS/polCoef.h"
@@ -44,6 +45,7 @@
 #include "AMPTOOLS_AMPS/DblRegge_FastPi.h"
 #include "AMPTOOLS_AMPS/omegapi_amplitude.h"
 #include "AMPTOOLS_AMPS/Vec_ps_refl.h"
+#include "AMPTOOLS_AMPS/Iso_ps_refl.h"
 #include "AMPTOOLS_AMPS/PhaseOffset.h"
 #include "AMPTOOLS_AMPS/ComplexCoeff.h"
 #include "AMPTOOLS_AMPS/OmegaDalitz.h"
@@ -456,6 +458,7 @@ int main( int argc, char* argv[] ){
    AmpToolsInterface::registerAmplitude( DblRegge_FastPi() );
    AmpToolsInterface::registerAmplitude( omegapi_amplitude() );
    AmpToolsInterface::registerAmplitude( Vec_ps_refl() );
+   AmpToolsInterface::registerAmplitude( Iso_ps_refl() );
    AmpToolsInterface::registerAmplitude( PhaseOffset() );
    AmpToolsInterface::registerAmplitude( ComplexCoeff() );
    AmpToolsInterface::registerAmplitude( OmegaDalitz() );
@@ -469,7 +472,8 @@ int main( int argc, char* argv[] ){
    AmpToolsInterface::registerAmplitude( KopfKMatrixRho() );
    AmpToolsInterface::registerAmplitude( KopfKMatrixPi1() );
    AmpToolsInterface::registerAmplitude( Vec_ps_moment() );
-
+   AmpToolsInterface::registerAmplitude( PiPiSWaveAMPK() );
+   
    AmpToolsInterface::registerDataReader( ROOTDataReader() );
    AmpToolsInterface::registerDataReader( ROOTDataReaderBootstrap() );
    AmpToolsInterface::registerDataReader( ROOTDataReaderWithTCut() );
diff --git a/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc b/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc
index 46883d57e..7c73f2f73 100644
--- a/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc
+++ b/src/programs/AmplitudeAnalysis/fitMPI/fitMPI.cc
@@ -32,6 +32,7 @@
 #include "AMPTOOLS_AMPS/VecRadiative_SDME.h"
 #include "AMPTOOLS_AMPS/Zlm.h"
 #include "AMPTOOLS_AMPS/BreitWigner.h"
+#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h"
 #include "AMPTOOLS_AMPS/BreitWigner3body.h"
 #include "AMPTOOLS_AMPS/b1piAngAmp.h"
 #include "AMPTOOLS_AMPS/Uniform.h"
@@ -40,6 +41,8 @@
 #include "AMPTOOLS_AMPS/DblRegge_FastPi.h"
 #include "AMPTOOLS_AMPS/omegapi_amplitude.h"
 #include "AMPTOOLS_AMPS/Vec_ps_refl.h"
+#include "AMPTOOLS_AMPS/VecPs_3pi_refl.h"
+#include "AMPTOOLS_AMPS/Iso_ps_refl.h"
 #include "AMPTOOLS_AMPS/Piecewise.h"
 #include "AMPTOOLS_AMPS/Flatte.h"
 #include "AMPTOOLS_AMPS/PhaseOffset.h"
@@ -66,7 +69,7 @@ using std::complex;
 using namespace std;
 
 int rank_mpi;
-int size;
+int mpi_size;
 
 double runSingleFit(ConfigurationInfo* cfgInfo, bool useMinos, bool hesse, int maxIter, string seedfile, int eMatrixRequirement) {
    AmpToolsInterfaceMPI ati( cfgInfo );
@@ -325,7 +328,7 @@ int main( int argc, char* argv[] ){
    MPI_Init( &argc, &argv );
 
    MPI_Comm_rank( MPI_COMM_WORLD, &rank_mpi );
-   MPI_Comm_size( MPI_COMM_WORLD, &size );
+   MPI_Comm_size( MPI_COMM_WORLD, &mpi_size );
 
    // set default parameters
 
@@ -428,6 +431,8 @@ int main( int argc, char* argv[] ){
    AmpToolsInterface::registerAmplitude( DblRegge_FastPi() );
    AmpToolsInterface::registerAmplitude( omegapi_amplitude() );
    AmpToolsInterface::registerAmplitude( Vec_ps_refl() );
+   AmpToolsInterface::registerAmplitude( Iso_ps_refl() );
+   AmpToolsInterface::registerAmplitude( VecPs_3pi_refl() );
    AmpToolsInterface::registerAmplitude( Piecewise() );
    AmpToolsInterface::registerAmplitude( Flatte() );
    AmpToolsInterface::registerAmplitude( PhaseOffset() );
@@ -442,7 +447,9 @@ int main( int argc, char* argv[] ){
    AmpToolsInterface::registerAmplitude( KopfKMatrixRho() );
    AmpToolsInterface::registerAmplitude( KopfKMatrixPi1() );
    AmpToolsInterface::registerAmplitude( Vec_ps_moment() );
+   AmpToolsInterface::registerAmplitude( PiPiSWaveAMPK() );
 
+   
    AmpToolsInterface::registerDataReader( DataReaderMPI<ROOTDataReader>() );
    AmpToolsInterface::registerDataReader( DataReaderMPI<ROOTDataReaderBootstrap>() );
    AmpToolsInterface::registerDataReader( DataReaderMPI<ROOTDataReaderWithTCut>() );
diff --git a/src/programs/AmplitudeAnalysis/isops_plotter/SConscript b/src/programs/AmplitudeAnalysis/isops_plotter/SConscript
new file mode 100644
index 000000000..7dcb10d85
--- /dev/null
+++ b/src/programs/AmplitudeAnalysis/isops_plotter/SConscript
@@ -0,0 +1,22 @@
+
+import os
+import sbms
+
+# get env object and clone it
+Import('*')
+
+# Verify AMPTOOLS environment variable is set
+if os.getenv('AMPTOOLS', 'nada')!='nada' and os.getenv('AMPPLOTTER', 'nada')!='nada':
+
+   env = env.Clone()
+
+   AMPTOOLS_LIBS = "AMPTOOLS_AMPS AMPTOOLS_DATAIO AMPTOOLS_MCGEN UTILITIES"
+   env.AppendUnique(LIBS = AMPTOOLS_LIBS.split())
+
+   sbms.AddHDDM(env)
+   sbms.AddAmpTools(env)
+   sbms.AddAmpPlotter(env)
+   sbms.AddUtilities(env)
+   sbms.AddROOT(env)
+
+   sbms.executable(env)
diff --git a/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc
new file mode 100644
index 000000000..2edc2b6d3
--- /dev/null
+++ b/src/programs/AmplitudeAnalysis/isops_plotter/isops_plotter.cc
@@ -0,0 +1,500 @@
+#include <iostream>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <map>
+
+#include "TClass.h"
+#include "TApplication.h"
+#include "TGClient.h"
+#include "TROOT.h"
+#include "TH1.h"
+#include "TStyle.h"
+#include "TClass.h"
+#include "TFile.h"
+#include "TMultiGraph.h"
+#include "TGraphErrors.h"
+
+#include "IUAmpTools/AmpToolsInterface.h"
+#include "IUAmpTools/FitResults.h"
+
+#include "AmpPlotter/PlotterMainWindow.h"
+#include "AmpPlotter/PlotFactory.h"
+
+#include "AMPTOOLS_DATAIO/IsoPsPlotGenerator.h"
+#include "AMPTOOLS_DATAIO/ROOTDataReader.h"
+#include "AMPTOOLS_DATAIO/ROOTDataReaderBootstrap.h"
+#include "AMPTOOLS_DATAIO/ROOTDataReaderTEM.h"
+#include "AMPTOOLS_DATAIO/FSRootDataReader.h"
+#include "AMPTOOLS_AMPS/BreitWigner.h"
+#include "AMPTOOLS_AMPS/Uniform.h"
+#include "AMPTOOLS_AMPS/PiPiSWaveAMPK.h"
+#include "AMPTOOLS_AMPS/Iso_ps_refl.h"
+#include "AMPTOOLS_AMPS/PhaseOffset.h"
+#include "AMPTOOLS_AMPS/ComplexCoeff.h"
+#include "AMPTOOLS_AMPS/Piecewise.h"
+#include "AMPTOOLS_AMPS/OmegaDalitz.h"
+
+
+
+#include "MinuitInterface/MinuitMinimizationManager.h"
+#include "IUAmpTools/ConfigFileParser.h"
+#include "IUAmpTools/ConfigurationInfo.h"
+
+
+
+void atiSetup(){
+  
+  AmpToolsInterface::registerAmplitude( BreitWigner() );
+  AmpToolsInterface::registerAmplitude( Uniform() );
+  AmpToolsInterface::registerAmplitude( Iso_ps_refl() );
+  AmpToolsInterface::registerAmplitude( PhaseOffset() );
+  AmpToolsInterface::registerAmplitude( ComplexCoeff() );
+  AmpToolsInterface::registerAmplitude( Piecewise() );
+  AmpToolsInterface::registerAmplitude( OmegaDalitz() );
+  AmpToolsInterface::registerAmplitude( PiPiSWaveAMPK() );
+
+  
+  AmpToolsInterface::registerDataReader( ROOTDataReader() );
+  AmpToolsInterface::registerDataReader( FSRootDataReader() );
+  AmpToolsInterface::registerDataReader( ROOTDataReaderTEM() );
+}
+
+using namespace std;
+
+int main( int argc, char* argv[] ){
+
+
+    // ************************
+    // usage
+    // ************************
+
+  cout << endl << " *** Viewing Results Using AmpPlotter and writing root histograms *** " << endl << endl;
+
+  if (argc < 2){
+    cout << "Usage:" << endl << endl;
+    cout << "\tisops_plotter <results file name> -o <output root file name> -t <output text file name>" << endl << endl;
+    return 0;
+  }
+
+  bool showGui = false;
+  string resultsName(argv[1]);
+  string outName = "isops_histos.root";
+  string outparsName = "isops__fitpars.txt";
+  
+  for (int i = 2; i < argc; i++){
+
+    string arg(argv[i]);
+
+    if (arg == "-o")  outName = argv[++i];
+    if (arg == "-t")  outparsName = argv[++i];
+    if (arg == "-g")  showGui = true;
+
+    
+    if (arg == "-h"){
+      cout << endl << " Usage for: " << argv[0] << endl << endl;
+      cout << "\t -o <file>\t output root file path" << endl;
+      cout << "\t -t <file>\t output text file path" << endl;
+      cout << "\t -g <file>\t show GUI" << endl;  
+      exit(1);
+    }
+  }
+
+
+    // ************************
+    // parse the command line parameters
+    // ************************
+
+  cout << "Fit results file name    = " << resultsName << endl;
+  cout << "Output root file name    = " << outName << endl;
+  cout << "Output text file name    = " << outparsName << endl << endl;
+
+    // ************************
+    // load the results and display the configuration info
+    // ************************
+
+  cout << "Loading Fit results" << endl;
+  FitResults results( resultsName );
+  if( !results.valid() ){
+    
+    cout << "Invalid fit results in file:  " << resultsName << endl;
+    exit( 1 );
+  }
+  cout << "Fit results loaded" << endl;
+    // ************************
+    // set up the plot generator
+    // ************************
+	cout << "before atisetup();"<< endl;
+	atiSetup();
+        cout << "Plotgen results"<< endl;
+
+	//	IsoPsPlotGenerator plotGen( results, PlotGenerator::kNoGenMC ); // optional can be omitted
+	IsoPsPlotGenerator plotGen( results ); 
+	cout << " Initialized ati and PlotGen" << endl;
+
+ 
+    // *************************
+    // Define Amplitudes and Coherent sums
+    // *****************************
+
+    vector<string> reflname = {"Uniform","PosRefl", "NegRefl"};
+    vector<string> amphistname = {"Flat"};    
+    vector<string> pipiIsobar_amps = {"pipiIso_0-S","pipiIso_1+P-","pipiIso_1+P0","pipiIso_1+P+"};
+    vector<string> rhoIsobar_amps = {"rhoIso_0-P", "rhoIso_1+S-", "rhoIso_1+S0", "rhoIso_1+S+", "rhoIso_1+D-", "rhoIso_1+D0", "rhoIso_1+D+", "rhoIso_2+D--", "rhoIso_2+D-", "rhoIso_2+D0", "rhoIso_2+D+", "rhoIso_2+D++"}; 
+    vector<string> f2Isobar_amps = {"f2Iso_1+P-","f2Iso_1+P0","f2Iso_1+P+","f2Iso_2+P--","f2Iso_2+P-","f2Iso_2+P0","f2Iso_2+P+","f2Iso_2+P++","f2Iso_2-S--","f2Iso_2-S-","f2Iso_2-S0","f2Iso_2-S+","f2Iso_2-S++","f2Iso_2-D--","f2Iso_2-D-","f2Iso_2-D0","f2Iso_2-D+","f2Iso_2-D++"};
+
+    amphistname.insert(amphistname.end(), pipiIsobar_amps.begin(), pipiIsobar_amps.end());
+    amphistname.insert(amphistname.end(), rhoIsobar_amps.begin(), rhoIsobar_amps.end());
+    amphistname.insert(amphistname.end(), f2Isobar_amps.begin(), f2Isobar_amps.end());
+
+    vector<string> ampsumname = {"pipiIso_0-S","pipiIso_1+P","rhoIso_0-P","rhoIso_1+S","rhoIso_1+D","rhoIso_2+D","f2Iso_1+P","f2Iso_2+P","f2Iso_2-S","f2Iso_2-D"};
+
+    
+    
+    // ************************
+    // set up an output ROOT file to store histograms
+    // ************************
+
+     TFile* plotfile = new TFile( outName.c_str(), "recreate");
+     TH1::AddDirectory(kFALSE);
+     TDirectory* plotfiledir = plotfile->mkdir("Contributions");  
+
+     
+     // *************************
+     // Loop over different polarization types 
+    // *************************
+     size_t nReactions = results.reactionList().size();
+     std::map<std::string, TH1*> summedHists;
+
+     
+  for (unsigned int polType = 0; polType < nReactions; polType++) {
+
+    
+    string reactionName = results.reactionList()[polType];
+
+    //WRITING SEPARATE FILES FOR DIFF POLS    
+    //  outName = reactionName + ".root";
+    //TFile* plotfile = new TFile( outName.c_str(), "recreate");
+    //TH1::AddDirectory(kFALSE);
+
+    //    string dirname = "Contributions_"+reactionName;
+    //    TDirectory* plotfiledir = plotfile->mkdir(dirname.c_str());  
+
+  
+    plotGen.enableReaction( reactionName );
+    vector<string> sums = plotGen.uniqueSums();
+    vector<string> amps = plotGen.uniqueAmplitudes();
+    cout << "Reaction " << reactionName << " enabled with " << sums.size() << " sums and " << amps.size() << " amplitudes" << endl;
+
+    string locampname;
+  
+    // loop over sum configurations (one for each of the individual contributions, and the combined sum of all)
+    for (unsigned int irefl = 0; irefl <= reflname.size(); irefl++){
+
+      // loop over desired amplitudes 
+      for (unsigned int iamp = 0; iamp <= amphistname.size(); iamp++ ) {
+
+	// turn on all ampltiudes by default 
+	for (unsigned int jamp = 0; jamp < amps.size(); jamp++ ) plotGen.enableAmp( jamp );
+      
+	if (iamp < amphistname.size()){
+	  // amplitude name
+	  locampname = amphistname[iamp];	
+
+	  // turn on all sums by default
+	  for (unsigned int jsum = 0; jsum < sums.size(); jsum++) plotGen.enableSum(jsum);
+	
+	  //Arrange the negative and positive reflectivity sums 
+	  //Negative reflectivity: "ImagNegSign" & "RealPosSign", as in the config file
+	  //Positive reflectivity: "RealNegSign" & "ImagPosSign", as in the config file
+	   	
+	  if (irefl < 3) {
+	    for (unsigned int i = 0; i < sums.size(); i++){
+	      if( reflname[irefl] == "NegRefl" ){
+		if(sums[i].find("RealNegSign") != std::string::npos || sums[i].find("ImagPosSign") != std::string::npos || sums[i].find("UniBG") != std::string::npos){
+		  plotGen.disableSum(i);
+		  //cout<<"disable sum "<<sums[i]<<"\n";
+		}
+	      }
+	      if( reflname[irefl] == "PosRefl" ){
+		if(sums[i].find("ImagNegSign") != std::string::npos || sums[i].find("RealPosSign") != std::string::npos || sums[i].find("UniBG") != std::string::npos) {
+		  //cout<<"disable sum "<<sums[i]<<"\n";
+		  plotGen.disableSum(i);
+		}
+	      }
+	      if( reflname[irefl] == "Uniform" ){
+		if(sums[i].find("ImagNegSign") != std::string::npos || sums[i].find("RealPosSign") != std::string::npos
+		   || sums[i].find("RealNegSign") != std::string::npos || sums[i].find("ImagPosSign") != std::string::npos) {
+		  //cout<<"disable sum "<<sums[i]<<"\n";
+		  plotGen.disableSum(i);
+		}
+	      }
+	    }	 
+	  }
+	  
+	  //Pick up the right amplitude by excluding all nonmatching amplitude names
+	  for (unsigned int jamp = 0; jamp < amps.size(); jamp++ ) 
+	    if( amps[jamp].find(locampname.data()) == std::string::npos )  plotGen.disableAmp( jamp );
+
+	} // close if (iamp < amphistname.size())
+
+
+	
+	// GENERATE HISTOGRAMS BY MEANS OF PLOT GENERATOR	
+    
+      cout << "Looping over input data" << endl;
+      // loop over data, accMC, and genMC
+      for (unsigned int iplot = 0; iplot < PlotGenerator::kNumTypes; iplot++){
+		bool singleData =  irefl == reflname.size() && iamp == amphistname.size();
+		bool singleFlatWave = (irefl == 0 && iamp > 0) || (irefl > 0 && iamp == 0);
+		
+		// if (iplot == PlotGenerator::kGenMC) continue; // no acceptance correction
+		if ( iplot == PlotGenerator::kData && !singleData ) continue; // only plot data once
+		if ( iplot == PlotGenerator::kBkgnd && !singleData ) continue; // only plot background once
+		if ( iplot == PlotGenerator::kAccMC && singleFlatWave ) continue; // only plot Flat wave once
+		if ( iplot == PlotGenerator::kGenMC && singleFlatWave ) continue; // only plot Flat wave once
+		
+		
+	// loop over different variables
+	for (unsigned int ivar  = 0; ivar  < IsoPsPlotGenerator::kNumHists; ivar++){
+	  
+	  // set reaction name as a prefix to the histogram name if all reactions are saved separately
+	  // set a common name or no name as a prefix to the histogram name if all reactions are summed up 
+	  //string histname = reactionName;
+	  string histname = "";
+	  
+	  
+	  if (ivar == IsoPsPlotGenerator::kProd_Ang)  histname += "Prod_Ang";
+	  else if (ivar == IsoPsPlotGenerator::kCosTheta)  histname += "CosTheta_GJ";
+	  else if (ivar == IsoPsPlotGenerator::kPhi)  histname += "Phi_GJ";
+	  else if (ivar == IsoPsPlotGenerator::kCosThetaH)  histname += "CosTheta_HF";
+	  else if (ivar == IsoPsPlotGenerator::kPhiH)  histname += "Phi_HF";
+	  else if (ivar == IsoPsPlotGenerator::kIsoMass)  histname += "MIso";
+	  else if (ivar == IsoPsPlotGenerator::kIsoPsMass)  histname += "MIsoPs";
+	  else if (ivar == IsoPsPlotGenerator::kt)  histname += "minust";
+	  else if (ivar == IsoPsPlotGenerator::kRecoilMass)  histname += "ProtonPiplusL_M";
+	  else if (ivar == IsoPsPlotGenerator::kProtonPsMass)  histname += "ProtonPiminus_M";
+	  else if (ivar == IsoPsPlotGenerator::kRecoilPsMass)  histname += "ProtonPiplusLPiminus_M";
+	  else continue;	  
+
+	  if (iplot == PlotGenerator::kData) histname += "_data";
+	  if (iplot == PlotGenerator::kBkgnd) histname += "_bkgnd";
+	  if (iplot == PlotGenerator::kAccMC) histname += "_fit";
+	  if (iplot == PlotGenerator::kGenMC) histname += "_gen";
+	  
+	  if (irefl < reflname.size()){
+	    // get name of sum for naming histogram
+	    string sumName = reflname[irefl];
+	    histname += "_";
+	    histname += sumName;
+	  }
+	  //	  if (iamp > 0 && iamp < amphistname.size()) {
+	  if (iamp < amphistname.size()) {
+	    // get name of amp for naming histogram  
+	    histname += "_";
+	    histname += amphistname[iamp];
+	  }
+	  
+	  Histogram* hist = plotGen.projection(ivar, reactionName, iplot);
+	  TH1* thist = hist->toRoot();
+	  thist->SetName(histname.c_str());
+
+	  // sum histograms across reactions
+	  if (summedHists.find(histname) == summedHists.end())
+	    summedHists[histname] = (TH1*) thist->Clone();
+	  else 
+	    summedHists[histname]->Add(thist);
+
+	  //write separate histogram for each reaction
+	  //if (iamp > 0 && iamp < amphistname.size())  plotfiledir->cd();
+	  //else  plotfile->cd();
+	  //thist->Write();
+	  
+	  delete thist;
+	}
+      }
+    } // end of loop over amplitudes
+   } // end of loop over sum configurations
+  }// end of loop over 'reactions'
+
+  
+
+  //Now, write (nonzero) summed up histograms to the root file
+  for (auto it = summedHists.begin(); it != summedHists.end(); it++) {
+
+    std::string hsummed_name = it->first;
+    TH1* hsummed = it->second;    
+
+    if (hsummed->Integral() == 0.) continue;
+    else hsummed->SetName(hsummed_name.c_str());
+    
+    if (hsummed_name.find("+") != std::string::npos || hsummed_name.find("-") != std::string::npos)  plotfiledir->cd();              
+    else  plotfile->cd();    
+    hsummed->Write();
+  }
+
+  
+  plotfile->Close(); 
+
+
+
+  // CALCULATE INTENSITY FRACTIONS AND PHASE DIFFERENCES
+  // All 'reactions' are already turned on
+  // The next calculations only need to happen one time
+  
+  // model parameters
+  // cout << "Checking Parameters" << endl;
+  
+  // parameters to check
+  vector< string > pars;  
+  //  pars.push_back("dsratio");
+  
+  //Text file for writing the parameters
+  ofstream outfile;
+  outfile.open( outparsName );
+
+  for(unsigned int i = 0; i<pars.size(); i++) {
+    double parValue = results.parValue( pars[i] );
+    double parError = results.parError( pars[i] );
+    outfile << parValue << "\t" << parError << "\t" << endl;
+  }
+
+
+  
+  //Define vector with full amplitudes
+  vector<string> fullamps = plotGen.fullAmplitudes();
+
+  //Define vectors for combining amplitudes with unique JPLMeps 
+  const int nAmps = amphistname.size();
+  vector<string> ampPosRefl[nAmps];
+  vector<string> ampNegRefl[nAmps];
+  vector< pair<string,string> > phaseDiffNames;
+
+  //Define vectors for combining amplitudes with unique JPLeps 
+  const int nSums = ampsumname.size();
+  vector<string> ampsumPosRefl[nSums];
+  vector<string> ampsumNegRefl[nSums];
+
+
+  for(unsigned int i = 0; i < fullamps.size(); i++){    
+
+    // Split by reflectivities and grab contributions for every JPLM amplitude 
+    for(int iamp=0; iamp<nAmps; iamp++) {
+	    string locampname = amphistname[iamp];
+	    
+	    if (fullamps[i].find("::" + locampname) == std::string::npos) continue;
+
+	    if (fullamps[i].find("ImagNegSign") != std::string::npos || fullamps[i].find("RealPosSign") != std::string::npos)
+	      ampNegRefl[iamp].push_back(fullamps[i]);
+
+	    if (fullamps[i].find("ImagPosSign") != std::string::npos || fullamps[i].find("RealNegSign") != std::string::npos)
+	      ampPosRefl[iamp].push_back(fullamps[i]);
+    }
+
+    // Split by reflectivities and grab contributions for every JPL sum 
+    for(int isum=0; isum<nSums; isum++) {
+	    string locsumname = ampsumname[isum];
+	    
+	    if (fullamps[i].find("::" + locsumname) == std::string::npos) continue;
+
+	    if (fullamps[i].find("ImagNegSign") != std::string::npos || fullamps[i].find("RealPosSign") != std::string::npos)
+	      ampsumNegRefl[isum].push_back(fullamps[i]);
+	      
+	    if (fullamps[i].find("ImagPosSign") != std::string::npos || fullamps[i].find("RealNegSign") != std::string::npos)
+	      ampsumPosRefl[isum].push_back(fullamps[i]);
+     }
+    
+    
+    // second loop over amplitudes to get phase difference names
+    for(unsigned int j = i+1; j < fullamps.size(); j++){
+
+      // leave only the Spring2017_PARA_0::ImagPosSign and Spring2017_PARA_0::ImagNegSign coherent sums
+      if (fullamps[i].find("Spring2017_PARA_0") == std::string::npos || fullamps[j].find("Spring2017_PARA_0") == std::string::npos) continue;      
+      if (fullamps[i].find("UniBG") != std::string::npos || fullamps[j].find("UniBG") != std::string::npos) continue;      
+      if (fullamps[i].find("Real") != std::string::npos || fullamps[j].find("Real") != std::string::npos) continue;      
+      
+      // only pair amplitudes in the same coherent sums
+      if (fullamps[i].find("ImagNegSign") != std::string::npos && fullamps[j].find("ImagNegSign") == std::string::npos) continue;
+      if (fullamps[i].find("ImagPosSign") != std::string::npos && fullamps[j].find("ImagPosSign") == std::string::npos) continue;
+	    
+      phaseDiffNames.push_back( std::make_pair(fullamps[i], fullamps[j]) );
+      //      std::cout << "PAIRED:\t" << fullamps[i] << "\t" << fullamps[j] << endl; 
+    }
+  }
+
+
+  outfile << "TOTAL EVENTS = " << results.intensity().first << " +- " << results.intensity().second << endl;
+  
+  cout<<"Computing intensities of amplitudes"<<endl;
+  for (unsigned int i = 0; i < fullamps.size(); i++){
+    vector<string> useamp;
+    useamp.push_back(fullamps[i]);
+    outfile << "FIT FRACTION " << fullamps[i] << " = " << results.intensity(useamp).first / results.intensity().first << " +- " << results.intensity(useamp).second / results.intensity().first <<  endl;
+  }
+  
+  cout<<"Computing phase differences"<<endl;
+  for(unsigned int i = 0; i < phaseDiffNames.size(); i++) {
+	  pair <double, double> phaseDiff = results.phaseDiff( phaseDiffNames[i].first, phaseDiffNames[i].second );
+	  outfile << "PHASE DIFF " << phaseDiffNames[i].first << " " << phaseDiffNames[i].second << " " << phaseDiff.first << " " << phaseDiff.second << endl;
+  }
+
+  cout<<"Computing intensities of coherent sums"<<endl;
+  for(int i = 0; i < nAmps; i++){
+    if(ampPosRefl[i].empty()) continue;
+    outfile << "FIT FRACTION (coherent sum) PosRefl " << amphistname[i] << " = " << results.intensity(ampPosRefl[i]).first / results.intensity().first << " +- "
+	    << results.intensity(ampPosRefl[i]).second / results.intensity().first << endl;
+     outfile << "FIT FRACTION (coherent sum) NegRefl " << amphistname[i] << " = " << results.intensity(ampNegRefl[i]).first / results.intensity().first << " +- "
+	     << results.intensity(ampNegRefl[i]).second / results.intensity().first << endl;
+  }
+
+  for(int i = 0; i < nSums; i++){
+    if(ampsumPosRefl[i].empty()) continue;
+    outfile << "FIT FRACTION (coherent sum) PosRefl " << ampsumname[i] << " = " << results.intensity(ampsumPosRefl[i]).first / results.intensity().first << " +- "
+	    << results.intensity(ampsumPosRefl[i]).second / results.intensity().first << endl;
+     outfile << "FIT FRACTION (coherent sum) NegRefl " << ampsumname[i] << " = " << results.intensity(ampsumNegRefl[i]).first / results.intensity().first << " +- "
+	     << results.intensity(ampsumNegRefl[i]).second / results.intensity().first << endl;
+  }
+
+
+
+  
+  
+  // covariance matrix
+  vector< vector< double > > covMatrix;
+  covMatrix = results.errorMatrix();
+
+
+
+    // ************************
+    // start the GUI
+    // ************************
+
+  if(showGui) {
+
+	  cout << ">> Plot generator ready, starting GUI..." << endl;
+	  
+	  int dummy_argc = 0;
+	  char* dummy_argv[] = {};  
+	  TApplication app( "app", &dummy_argc, dummy_argv );
+	  
+	  gStyle->SetFillColor(10);
+	  gStyle->SetCanvasColor(10);
+	  gStyle->SetPadColor(10);
+	  gStyle->SetFillStyle(1001);
+	  gStyle->SetPalette(1);
+	  gStyle->SetFrameFillColor(10);
+	  gStyle->SetFrameFillStyle(1001);
+   
+     cout << " Initialized App " << endl;     
+	  PlotFactory factory( plotGen );	
+     cout << " Created Plot Factory " << endl;
+	  PlotterMainWindow mainFrame( gClient->GetRoot(), factory );
+     cout << " Main frame created " << endl;
+     
+	  app.Run();
+     cout << " App running" << endl;
+  }
+    
+  return 0;
+
+}
diff --git a/src/programs/Simulation/gen_vec_ps/gen_vec_ps.cc b/src/programs/Simulation/gen_vec_ps/gen_vec_ps.cc
index 1205ca433..5e55f3fa6 100644
--- a/src/programs/Simulation/gen_vec_ps/gen_vec_ps.cc
+++ b/src/programs/Simulation/gen_vec_ps/gen_vec_ps.cc
@@ -567,7 +567,7 @@ int main( int argc, char* argv[] ){
 
 					TLorentzVector Gammap = beam + target;
 					if( polAngle == -1 ){
-						fprintf(stderr, "Error: do not use amorphous value (got %d)\n", polAngle);
+						fprintf(stderr, "Error: do not use amorphous value (got %f)\n", polAngle);
         				return -1; 
 					}
                     vector <double> xDecayAngles = getXDecayAngles(polAngle, beam, Gammap, isobar, resonance);