diff --git a/inc/LauGounarisSakuraiRes.hh b/inc/LauGounarisSakuraiRes.hh
index cefa227..fed2168 100644
--- a/inc/LauGounarisSakuraiRes.hh
+++ b/inc/LauGounarisSakuraiRes.hh
@@ -1,126 +1,140 @@
 
 /*
 Copyright 2006 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauGounarisSakuraiRes.hh
     \brief File containing declaration of LauGounarisSakuraiRes class.
 */
 
 /*! \class LauGounarisSakuraiRes
     \brief Class for defininf the Gounaris-Sakurai resonance model
 
     Class for defining the relativistic Gounaris-Sakurai resonance model, which
     includes the use of Blatt-Weisskopf barrier factors.
 */
 
 #ifndef LAU_GOUNARIS_SAKURAI_RES
 #define LAU_GOUNARIS_SAKURAI_RES
 
 #include "TString.h"
 
 #include "LauComplex.hh"
 #include "LauAbsResonance.hh"
 
 
 class LauGounarisSakuraiRes : public LauAbsResonance {
 
 	public:
 		//! Constructor
 		/*!
 			\param [in] resInfo the object containing information on the resonance name, mass, width, spin, charge, etc.
 			\param [in] resPairAmpInt the number of the daughter not produced by the resonance 
 			\param [in] daughters the daughter particles
 		*/	
 		LauGounarisSakuraiRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters);
 
 		//! Destructor
 		virtual ~LauGounarisSakuraiRes();
 
 		//! Initialise the model
 		virtual void initialise();
 
 		//! Get the resonance model type
                 /*!
                         \return the resonance model type
                 */
 		virtual LauAbsResonance::LauResonanceModel getResonanceModel() const {return LauAbsResonance::GS;}
 
 		//! Retrieve the resonance parameters, e.g. so that they can be loaded into a fit
 		/*!
 		    \return floating parameters of the resonance
 		*/
 		virtual const std::vector<LauParameter*>& getFloatingParameters();
 
+		//! Optionally disable scaling Gamma(m) by the Blatt-Weisskopf factor
+		/*!
+			Equation 11 of the original paper does not contain such
+			a factor, so we allow here the option to switch it off
+			to better match the original formulation.
+			It is switched on by default for backward compatibility.
+
+			\param [in] flag if true, disable the scaling, if false, revert to default behaviour
+		*/
+		static void disableBWScalingOfWidth( const Bool_t flag = kTRUE ) { disableBWScalingOfWidth_ = flag; }
+
 	protected:
 		//! Complex resonant amplitude
 		/*!
 			\param [in] mass appropriate invariant mass for the resonance
 			\param [in] spinTerm spin term
 		*/	
 		virtual LauComplex resAmp(Double_t mass, Double_t spinTerm);
 
 	private:
 		//! Copy constructor (not implemented)
 		LauGounarisSakuraiRes(const LauGounarisSakuraiRes& rhs);
 
 		//! Copy assignment operator (not implemented)
 		LauGounarisSakuraiRes& operator=(const LauGounarisSakuraiRes& rhs);
 
 		//! Momentum of the daughters in the resonance rest frame (at pole mass)
 		Double_t q0_; 
 		//! The resonance mass
 		Double_t resMass_;
 		//! Square of the resonance mass
 		Double_t resMassSq_;
 		//! The resonance width
 		Double_t resWidth_;
 		//! The resonance barrier radius
 		Double_t resRadius_;
 		//! The parent barrier radius
 		Double_t parRadius_;
 		//! Sum of the two daughter masses 
 		Double_t mDaugSum_; 
 		//! Square of the sum of the two daughter masses
 		Double_t mDaugSumSq_;
 		//! Difference between the two daughter masses
 		Double_t mDaugDiff_; 
 		//! Square of the difference of the two daughter masses
 		Double_t mDaugDiffSq_;
 		//! Square of the parent mass
 		Double_t mParentSq_; 
 		//! Square of the bachelor mass
 		Double_t mBachSq_;
 		//! Extra parameter required by GS shape
 		Double_t h0_; 
 		//! Extra parameter required by GS shape
 		Double_t dhdm0_; 
 		//! Extra parameter required by GS shape
 		Double_t d_;
 		//! Value of the form factor for resonance decay (at pole mass)
 		Double_t FR0_; 
 
+		//! Optionally disable scaling Gamma(m) by Blatt-Weisskopf factor
+		static Bool_t disableBWScalingOfWidth_;
+
 		ClassDef(LauGounarisSakuraiRes,0) // Gounaris-Sakurai resonance model
 
 };
 
 #endif
diff --git a/src/LauGounarisSakuraiRes.cc b/src/LauGounarisSakuraiRes.cc
index 1b5d780..3471e65 100644
--- a/src/LauGounarisSakuraiRes.cc
+++ b/src/LauGounarisSakuraiRes.cc
@@ -1,253 +1,254 @@
 
 /*
 Copyright 2006 University of Warwick
 
 Licensed under the Apache License, Version 2.0 (the "License");
 you may not use this file except in compliance with the License.
 You may obtain a copy of the License at
 
     http://www.apache.org/licenses/LICENSE-2.0
 
 Unless required by applicable law or agreed to in writing, software
 distributed under the License is distributed on an "AS IS" BASIS,
 WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 See the License for the specific language governing permissions and
 limitations under the License.
 */
 
 /*
 Laura++ package authors:
 John Back
 Paul Harrison
 Thomas Latham
 */
 
 /*! \file LauGounarisSakuraiRes.cc
     \brief File containing implementation of LauGounarisSakuraiRes class.
 */
 
 #include <iostream>
 
 #include "LauConstants.hh"
 #include "LauGounarisSakuraiRes.hh"
 
 ClassImp(LauGounarisSakuraiRes)
 
+Bool_t LauGounarisSakuraiRes::disableBWScalingOfWidth_{kFALSE};
 
 LauGounarisSakuraiRes::LauGounarisSakuraiRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
 	LauAbsResonance(resInfo, resPairAmpInt, daughters),
 	q0_(0.0),
 	resMass_(0.0),
 	resMassSq_(0.0),
 	resWidth_(0.0),
 	resRadius_(0.0),
 	parRadius_(0.0),
 	mDaugSum_(0.0),
 	mDaugSumSq_(0.0),
 	mDaugDiff_(0.0),
 	mDaugDiffSq_(0.0),
 	mParentSq_(0.0),
 	mBachSq_(0.0),
 	h0_(0.0),
 	dhdm0_(0.0),
 	d_(0.0),
 	FR0_(1.0)
 {
 }
 
 LauGounarisSakuraiRes::~LauGounarisSakuraiRes()
 {
 }
 
 void LauGounarisSakuraiRes::initialise()
 {
 	// Set-up various constants. This must be called again if the mass/width/spin
 	// of a resonance changes...  
 
 	resMass_ = this->getMass();
 	resWidth_ = this->getWidth();
 	resRadius_ = this->getResRadius();
 	parRadius_ = this->getParRadius();
 
 	Double_t massDaug1 = this->getMassDaug1();
 	Double_t massDaug2 = this->getMassDaug2();
 	Double_t massBachelor = this->getMassBachelor();
 	Double_t massParent = this->getMassParent();
 
 	// Check that the spin is 1
 	Int_t resSpin = this->getSpin();
 	if (resSpin != 1) {
 		std::cerr << "WARNING in LauGounarisSakuraiRes::initialise : Resonance spin is != 1. This lineshape is for the rho(770), setting the spin to 1." << std::endl;
 		this->changeResonance( -1.0, -1.0, 1 );
 		resSpin = this->getSpin();
 	}
 
 	// Create the mass squares, sums, differences etc.
 	resMassSq_ = resMass_*resMass_;
 	mDaugSum_  = massDaug1 + massDaug2;
 	mDaugSumSq_ = mDaugSum_*mDaugSum_;
 	mDaugDiff_ = massDaug1 - massDaug2;
 	mDaugDiffSq_ = mDaugDiff_*mDaugDiff_;
 	mParentSq_ = massParent*massParent;
 	mBachSq_ = massBachelor*massBachelor;
 
 	// Create an effective resonance pole mass to protect against resonances
 	// that are below threshold
 	Double_t effResMass = resMass_;
 	Double_t effResMassSq = resMassSq_;
 	if ( resMassSq_ - mDaugSumSq_ < 0.0 ) {
 		Double_t minMass = mDaugSum_;
 		Double_t maxMass = massParent - massBachelor;
 		Double_t tanhTerm = std::tanh( (resMass_ - ((minMass + maxMass)/2))/(maxMass-minMass));
 		effResMass = minMass + (maxMass-minMass)*(1+tanhTerm)/2;
 		effResMassSq = effResMass*effResMass;
 	}
 
 	// Decay momentum of either daughter in the resonance rest frame
 	// when resonance mass = rest-mass value, m_0 (PDG value)
 	Double_t term1 = effResMassSq - mDaugSumSq_;
 	Double_t term2 = effResMassSq - mDaugDiffSq_;
 	Double_t term12 = term1*term2;
 	if (term12 > 0.0) {
 		q0_ = TMath::Sqrt(term12)/(2.0*effResMass);
 	} else {
 		q0_ = 0.0;
 	}
 
 	// Calculate the Blatt-Weisskopf form factor for the case when m = m_0
 	FR0_ = 1.0;
 	if ( this->getSpin() > 0 ) {
 		const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
 		if ( resBWFactor != nullptr ) {
 			FR0_ = resBWFactor->calcFormFactor(q0_);
 		}
 	}
 
 	// Calculate the extra things needed by the G-S shape
 	h0_ = 2.0*LauConstants::invPi * q0_/effResMass * TMath::Log((effResMass + 2.0*q0_)/(2.0*LauConstants::mPi));
 	dhdm0_ = h0_ * (1.0/(8.0*q0_*q0_) - 1.0/(2.0*effResMassSq)) + 1.0/(LauConstants::twoPi*effResMassSq);
 	d_ = 3.0*LauConstants::invPi * LauConstants::mPi*LauConstants::mPi/(q0_*q0_)
 		* TMath::Log((effResMass + 2.0*q0_)/(2.0*LauConstants::mPi))
 		+ effResMass/(LauConstants::twoPi*q0_)
 		- LauConstants::mPi*LauConstants::mPi*effResMass/(LauConstants::pi*q0_*q0_*q0_);
 }
 
 LauComplex LauGounarisSakuraiRes::resAmp(Double_t mass, Double_t spinTerm)
 {
 	// This function returns the complex dynamical amplitude for a Breit-Wigner resonance,
 	// given the invariant mass and cos(helicity) values.
 
 	LauComplex resAmplitude(0.0, 0.0);
 
 	if (mass < 1e-10) {
 		std::cerr << "WARNING in LauGounarisSakuraiRes::amplitude : mass < 1e-10." << std::endl;
 		return LauComplex(0.0, 0.0);
 	} else if (q0_ < 1e-30) {
 		return LauComplex(0.0, 0.0);
 	}
 
 	// Calculate the width of the resonance (as a function of mass)
 	// First, calculate the various form factors.
 	// NB
 	// q is the momentum of either daughter in the resonance rest-frame,
 	// p is the momentum of the bachelor in the resonance rest-frame,
 	// pstar is the momentum of the bachelor in the parent rest-frame.
 	// These quantities have been calculate in LauAbsResonance::amplitude(...)
 
 	const Double_t resMass = this->getMass();
 	const Double_t resWidth = this->getWidth();
 	const Double_t resRadius = this->getResRadius();
 	const Double_t parRadius = this->getParRadius();
 
 	// If the mass is floating and its value has changed we need to
 	// recalculate everything that assumes that value
 	// Similarly for the BW radii
 	if ( ( (!this->fixMass()) && resMass != resMass_ ) ||
 	     ( (!this->fixResRadius()) && resRadius != resRadius_ ) ||
 	     ( (!this->fixParRadius()) && parRadius != parRadius_ ) ) {
 		this->initialise();
 	}
 
 	const Int_t resSpin = this->getSpin();
 	const Double_t q = this->getQ();
 	const Double_t p = this->getP();
 	const Double_t pstar = this->getPstar();
 
 	Double_t fFactorR(1.0);
 	Double_t fFactorB(1.0);
 	if ( resSpin > 0 ) {
 		const LauBlattWeisskopfFactor* resBWFactor = this->getResBWFactor();
 		if ( resBWFactor != nullptr ) {
 			fFactorR = resBWFactor->calcFormFactor(q);
 		}
 
 		const LauBlattWeisskopfFactor* parBWFactor = this->getParBWFactor();
 		if ( parBWFactor != nullptr ) {
 			switch ( parBWFactor->getRestFrame() ) {
 				case LauBlattWeisskopfFactor::ResonanceFrame:
 					fFactorB = parBWFactor->calcFormFactor(p);
 					break;
 				case LauBlattWeisskopfFactor::ParentFrame:
 					fFactorB = parBWFactor->calcFormFactor(pstar);
 					break;
 				case LauBlattWeisskopfFactor::Covariant:
 				{
 					Double_t covFactor = this->getCovFactor();
 					if ( resSpin > 2 ) {
 						covFactor = TMath::Power( covFactor, 1.0/resSpin );
 					} else if ( resSpin == 2 ) {
 						covFactor = TMath::Sqrt( covFactor );
 					}
 					fFactorB = parBWFactor->calcFormFactor(pstar*covFactor);
 					break;
 				}
 			}
 		}
 	}
 	const Double_t fFactorRRatio = fFactorR/FR0_;
 
 	const Double_t qRatio = q/q0_;
 	const Double_t qTerm = qRatio*qRatio*qRatio;
 
-	const Double_t totWidth = resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
+	const Double_t totWidth = disableBWScalingOfWidth_ ? resWidth*qTerm*(resMass/mass) : resWidth*qTerm*(resMass/mass)*fFactorRRatio*fFactorRRatio;
 
 	const Double_t massSq = mass*mass;
 	const Double_t massSqTerm = resMassSq_ - massSq;
 
 	const Double_t h = 2.0*LauConstants::invPi * q/mass * TMath::Log((mass + 2.0*q)/(2.0*LauConstants::mPi));
 	const Double_t f = resWidth * resMassSq_/(q0_*q0_*q0_) * (q*q * (h - h0_) + massSqTerm * q0_*q0_ * dhdm0_);
 
 	// Compute the complex amplitude
 	resAmplitude = LauComplex(massSqTerm + f, resMass*totWidth);
 
 	// Scale by the denominator factor, as well as the spin term and Blatt-Weisskopf factors
 	Double_t numerFactor = spinTerm*(1.0 + d_ * resWidth/resMass);
 	if (!this->ignoreBarrierScaling()) {
 		numerFactor *= fFactorR * fFactorB;
 	}
 	const Double_t denomFactor = (massSqTerm + f)*(massSqTerm + f) + resMassSq_*totWidth*totWidth;
 	resAmplitude.rescale(numerFactor/denomFactor);
 
 	return resAmplitude;
 }
 
 const std::vector<LauParameter*>& LauGounarisSakuraiRes::getFloatingParameters()
 {
 	this->clearFloatingParameters();
 
 	if ( ! this->fixMass() ) {
 		this->addFloatingParameter( this->getMassPar() );
 	}
 	if ( ! this->fixWidth() ) {
 		this->addFloatingParameter( this->getWidthPar() );
 	}
 	if ( ! this->fixResRadius() ) {
 		this->addFloatingParameter( this->getResBWFactor()->getRadiusParameter() );
 	}
 	if ( ! this->fixParRadius() ) {
 		this->addFloatingParameter( this->getParBWFactor()->getRadiusParameter() );
 	}
 
 	return this->getParameters();
 }