No OneTemporary
Actions

Size

83 KB

Subscribers

None

View Options

	diff --git a/examples/GenFit3K.cc b/examples/GenFit3K.cc
	index 130abf1..ffc8d0d 100644
	--- a/examples/GenFit3K.cc
	+++ b/examples/GenFit3K.cc
	@@ -1,191 +1,194 @@

	#include <cstdlib>
	#include <iostream>
	#include <vector>

	#include "TFile.h"
	#include "TH2.h"
	#include "TString.h"
	#include "TTree.h"

	#include "LauSimpleFitModel.hh"
	#include "LauBkgndDPModel.hh"
	#include "LauDaughters.hh"
	#include "LauEffModel.hh"
	#include "LauIsobarDynamics.hh"
	#include "LauMagPhaseCoeffSet.hh"
	#include "LauVetoes.hh"

	void usage( std::ostream& out, const TString& progName )
	{
	out<<"Usage:\n";
	out<<progName<<" gen [nExpt = 1] [firstExpt = 0]\n";
	out<<"or\n";
	out<<progName<<" fit <iFit> [nExpt = 1] [firstExpt = 0]"<<std::endl;
	}

	int main( int argc, char** argv )
	{
	// Process command-line arguments
	// Usage:
	// ./GenFit3pi gen [nExpt = 1] [firstExpt = 0]
	// or
	// ./GenFit3pi fit <iFit> [nExpt = 1] [firstExpt = 0]
	if ( argc < 2 ) {
	usage( std::cerr, argv[0] );
	return EXIT_FAILURE;
	}

	TString command = argv[1];
	command.ToLower();
	Int_t iFit(0);
	Int_t nExpt(1);
	Int_t firstExpt(0);
	if ( command == "gen" ) {
	if ( argc > 2 ) {
	nExpt = atoi( argv[2] );
	if ( argc > 3 ) {
	firstExpt = atoi( argv[3] );
	}
	}
	} else if ( command == "fit" ) {
	if ( argc < 3 ) {
	usage( std::cerr, argv[0] );
	return EXIT_FAILURE;
	}
	iFit = atoi( argv[2] );
	if ( argc > 3 ) {
	nExpt = atoi( argv[3] );
	if ( argc > 4 ) {
	firstExpt = atoi( argv[4] );
	}
	}
	} else {
	usage( std::cerr, argv[0] );
	return EXIT_FAILURE;
	}

	// If you want to use square DP histograms for efficiency,
	// backgrounds or you just want the square DP co-ordinates
	// stored in the toy MC ntuple then set this to kTRUE
	Bool_t squareDP = kFALSE;

	// This defines the DP => decay is B+ -> K+ K+ K-
	// Particle 1 = K+
	// Particle 2 = K+
	// Particle 3 = K-
	// The DP is defined in terms of m13Sq and m23Sq
	LauDaughters* daughters = new LauDaughters("B+", "K+", "K+", "K-", squareDP);

	// Optionally apply some vetoes to the DP
	LauVetoes* vetoes = new LauVetoes();

	// Define the efficiency model (defaults to unity everywhere)
	// Can optionally provide a histogram to model variation over DP
	// (example syntax given in commented-out section)
	LauEffModel* effModel = new LauEffModel(daughters, vetoes);
	//TFile *effHistFile = TFile::Open("histoFiles/effHistos.root", "read");
	//TH2* effHist = dynamic_cast<TH2*>(effHistFile->Get("effHist"));
	//Bool_t useInterpolation = kTRUE;
	//Bool_t fluctuateBins = kFALSE;
	//Bool_t useUpperHalf = kTRUE;
	//effModel->setEffHisto(effHist, useInterpolation, fluctuateBins, 0.0, 0.0, useUpperHalf, squareDP);

	// Create the isobar model
	LauIsobarDynamics* sigModel = new LauIsobarDynamics(daughters, effModel);

	// Add various components to the isobar model,
	// optionally allowing the masses and width to float in the fit
	LauAbsResonance* res(0);

	//addResonance arguments: resName, resPairAmpInt, resType
	res = sigModel->addResonance("phi(1020)", 1, LauAbsResonance::RelBW);
	//changeResonance arguments: newMass, newWidth, newSpin
	res->changeResonance(1.019, 0.0044, -1);
	res->fixMass(kFALSE);
	res->fixWidth(kFALSE);

	res = sigModel->addResonance("f'_2(1525)", 1, LauAbsResonance::RelBW);
	res->fixMass(kFALSE);
	res->fixWidth(kFALSE);

	res = sigModel->addResonance("NonReson", 0, LauAbsResonance::FlatNR);

	// Reset the maximum signal DP ASq value
	// This will be automatically adjusted to avoid bias or extreme
	// inefficiency if you get the value wrong but best to set this by
	// hand once you've found the right value through some trial and
	// error.
	sigModel->setASqMaxValue(16.4);

	// Create the fit model
	LauSimpleFitModel* fitModel = new LauSimpleFitModel(sigModel);

	// Create the complex coefficients for the isobar model
	// Here we're using the magnitude and phase form:
	// c_j = a_j exp(i*delta_j)
	std::vector<LauAbsCoeffSet*> coeffset;
	coeffset.push_back( new LauMagPhaseCoeffSet("phi(1020)", 1.0, 0.0, kTRUE, kTRUE) );
	coeffset.push_back( new LauMagPhaseCoeffSet("f'_2(1525)", 1.0, 0.0, kFALSE, kFALSE) );
	coeffset.push_back( new LauMagPhaseCoeffSet("NonReson", 1.0, 0.0, kFALSE, kFALSE) );

	for (std::vector<LauAbsCoeffSet*>::iterator iter=coeffset.begin(); iter!=coeffset.end(); ++iter) {
	fitModel->setAmpCoeffSet(*iter);
	}

	// Set the signal yield and define whether it is fixed or floated
	Int_t nSigEvents = 5000;
	Bool_t fixNSigEvents = kFALSE;
	LauParameter * signalEvents = new LauParameter("signalEvents", nSigEvents, -1.0nSigEvents, 2.0nSigEvents, fixNSigEvents);
	fitModel->setNSigEvents(signalEvents);

	// Set the number of experiments to generate or fit and which
	// experiment to start with
	fitModel->setNExpts( nExpt, firstExpt );

	// Switch on/off calculation of asymmetric errors.
	fitModel->useAsymmFitErrors(kFALSE);

	// Randomise initial fit values for the signal mode
	fitModel->useRandomInitFitPars(kFALSE);

	// Switch on/off Poissonian smearing of total number of events
	fitModel->doPoissonSmearing(kTRUE);

	// Switch on/off Extended ML Fit option
	Bool_t emlFit = ( fitModel->nBkgndClasses() > 0 );
	fitModel->doEMLFit(emlFit);

	+ // Switch on the two-stage fit (for the resonance parameters)
	+ fitModel->twoStageFit(kTRUE);
	+
	TString dataFile("data.root");

	TString treeName("genResults");
	TString rootFileName("");
	TString tableFileName("");
	TString fitToyFileName("fitToyMC_");
	TString splotFileName("splot_");
	if (command == "fit") {
	rootFileName = "fit"; rootFileName += iFit;
	rootFileName += "_expt_"; rootFileName += firstExpt;
	rootFileName += "-"; rootFileName += (firstExpt+nExpt-1);
	rootFileName += ".root";
	tableFileName = "fitResults_"; tableFileName += iFit;
	fitToyFileName += iFit;
	fitToyFileName += ".root";
	splotFileName += iFit;
	splotFileName += ".root";
	} else {
	rootFileName = "dummy.root";
	tableFileName = "genResults";
	}

	// Generate toy from the fitted parameters
	//fitModel->compareFitData(10, fitToyFileName);

	// Write out per-event likelihoods and sWeights
	//fitModel->writeSPlotData(splotFileName, "splot", kFALSE);

	// Execute the generation/fit
	fitModel->run( command, dataFile, treeName, rootFileName, tableFileName );

	return EXIT_SUCCESS;
	}
	diff --git a/src/LauBelleNR.cc b/src/LauBelleNR.cc
	index 15284be..37a0bfd 100644
	--- a/src/LauBelleNR.cc
	+++ b/src/LauBelleNR.cc
	@@ -1,132 +1,133 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauBelleNR.cc
	\brief File containing implementation of LauBelleNR class.
	*/

	#include <iostream>

	#include "TMath.h"

	#include "LauBelleNR.hh"
	#include "LauDaughters.hh"
	#include "LauParameter.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauBelleNR)


	LauBelleNR::LauBelleNR(LauResonanceInfo* resInfo, const LauAbsResonance::LauResonanceModel resType,
	const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	alpha_(0),
	model_(resType)
	{
	TString parName = this->getSanitisedName();
	parName += "_alpha";
	alpha_ = resInfo->getExtraParameter( parName );
	if ( alpha_ == 0 ) {
	alpha_ = new LauParameter( parName, 0.0, 0.0, 10.0, kTRUE );
	+ alpha_->secondStage(kTRUE);
	resInfo->addExtraParameter( alpha_ );
	}
	}

	LauBelleNR::~LauBelleNR()
	{
	delete alpha_;
	}

	void LauBelleNR::initialise()
	{
	const LauDaughters* daughters = this->getDaughters();
	Int_t resPairAmpInt = this->getPairInt();
	if ( daughters->gotSymmetricalDP() && resPairAmpInt != 3 ) {
	std::cerr << "WARNING in LauBelleNR::initialise : Dalitz plot is symmetric - this lineshape is not appropriate." << std::endl;
	}

	if ( model_ != LauAbsResonance::BelleNR && model_ != LauAbsResonance::PowerLawNR ) {
	std::cerr << "WARNING in LauBelleNR::initialise : Unknown model requested, defaulting to exponential." << std::endl;
	model_ = LauAbsResonance::BelleNR;
	}

	// TODO - should the BelleNR ignore the momenta in the spin term?
	}

	LauComplex LauBelleNR::resAmp(Double_t mass, Double_t spinTerm)
	{
	Double_t magnitude(1.0);

	Double_t alpha = this->getAlpha();

	if ( model_ == LauAbsResonance::BelleNR ) {
	magnitude = spinTerm * TMath::Exp(-alphamassmass);
	} else if ( model_ == LauAbsResonance::PowerLawNR ) {
	magnitude = spinTerm * TMath::Power(mass*mass, -alpha);
	}

	LauComplex resAmplitude(magnitude, 0.0);

	return resAmplitude;
	}

	const std::vector<LauParameter*>& LauBelleNR::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixAlpha() ) {
	this->addFloatingParameter( alpha_ );
	}

	return this->getParameters();
	}

	void LauBelleNR::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the lineshape
	if (name == "alpha") {
	this->setAlpha(value);
	std::cout << "INFO in LauBelleNR::setResonanceParameter : Setting parameter alpha = " << this->getAlpha() << std::endl;
	} else {
	std::cerr << "WARNING in LauBelleNR::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauBelleNR::floatResonanceParameter(const TString& name)
	{
	if (name == "alpha") {
	if ( alpha_->fixed() ) {
	alpha_->fixed( kFALSE );
	this->addFloatingParameter( alpha_ );
	} else {
	std::cerr << "WARNING in LauBelleNR::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauBelleNR::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauBelleNR::getResonanceParameter(const TString& name)
	{
	if (name == "alpha") {
	return alpha_;
	} else {
	std::cerr << "WARNING in LauBelleNR::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauBelleNR::setAlpha(const Double_t alpha)
	{
	alpha_->value( alpha );
	alpha_->genValue( alpha );
	alpha_->initValue( alpha );
	}

	diff --git a/src/LauBelleSymNR.cc b/src/LauBelleSymNR.cc
	index f7c1ee7..2227171 100644
	--- a/src/LauBelleSymNR.cc
	+++ b/src/LauBelleSymNR.cc
	@@ -1,156 +1,157 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauBelleSymNR.cc
	\brief File containing implementation of LauBelleSymNR class.
	*/

	#include <iostream>

	#include "TMath.h"

	#include "LauBelleSymNR.hh"
	#include "LauDaughters.hh"
	#include "LauKinematics.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauBelleSymNR)


	LauBelleSymNR::LauBelleSymNR(LauResonanceInfo* resInfo, const LauAbsResonance::LauResonanceModel resType,
	const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	alpha_(0),
	model_(resType)
	{
	TString parName = this->getSanitisedName();
	parName += "_alpha";
	alpha_ = resInfo->getExtraParameter( parName );
	if ( alpha_ == 0 ) {
	alpha_ = new LauParameter( parName, 0.0, 0.0, 10.0, kTRUE );
	+ alpha_->secondStage(kTRUE);
	resInfo->addExtraParameter( alpha_ );
	}
	}

	LauBelleSymNR::~LauBelleSymNR()
	{
	delete alpha_;
	}

	void LauBelleSymNR::initialise()
	{
	const LauDaughters* daughters = this->getDaughters();
	Int_t resPairAmpInt = this->getPairInt();
	if ( ! daughters->gotSymmetricalDP() ) {
	std::cerr << "WARNING in LauBelleSymNR::initialise : Dalitz plot is symmetric - this lineshape is not appropriate." << std::endl;
	}
	if ( resPairAmpInt == 3 ) {
	std::cerr << "WARNING in LauBelleSymNR::initialise : This lineshape is intended to be on the symmetrised axes of the DP." << std::endl;
	}

	if ( model_ != LauAbsResonance::BelleSymNR && model_ != LauAbsResonance::TaylorNR ) {
	std::cerr << "WARNING in LauBelleSymNR::initialise : Unknown model requested, defaulting to exponential." << std::endl;
	model_ = LauAbsResonance::BelleSymNR;
	}
	}

	LauComplex LauBelleSymNR::amplitude(const LauKinematics* kinematics)
	{
	// This function returns the complex dynamical amplitude for a Belle Non-Resonant distribution

	// Calculate for symmetric DPs, e.g. 3pi or 3K, by using shapeNo = 1 or 2
	// Have s<->t symmetry already done in Dynamics flip function.
	// For Kpipi or similar plots, one can use the separate exponentials
	// and consider them as two separate components with their own mag and phase.
	// For this shapeNo = 3 and shapeNo = 4 need to be used to create the two
	// individual amplitudes (with the same value of alpha).

	// Calculate Mandelstam variables.
	// s = m_13^2, t = m_23^2, u = m_12^2.
	Double_t s = kinematics->getm13Sq();
	Double_t t = kinematics->getm23Sq();
	//Double_t u = kinematics->getm12Sq();

	Double_t magnitude(1.0);

	Double_t alpha = this->getAlpha();

	if ( model_ == LauAbsResonance::BelleSymNR ) {
	magnitude = TMath::Exp(-alphas) + TMath::Exp(-alphat);
	} else if ( model_ == LauAbsResonance::TaylorNR ) {
	Double_t mParentSq = kinematics->getmParentSq();
	magnitude = alpha*(s + t)/mParentSq + 1.0;
	}

	LauComplex resAmplitude(magnitude, 0.0);

	return resAmplitude;
	}

	LauComplex LauBelleSymNR::resAmp(Double_t mass, Double_t spinTerm)
	{
	std::cerr << "ERROR in LauBelleSymNR : This method should never be called." << std::endl;
	std::cerr << " : Returning zero amplitude for mass = " << mass << " and spinTerm = " << spinTerm << "." << std::endl;
	return LauComplex(0.0, 0.0);
	}

	const std::vector<LauParameter*>& LauBelleSymNR::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixAlpha() ) {
	this->addFloatingParameter( alpha_ );
	}

	return this->getParameters();
	}

	void LauBelleSymNR::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the lineshape
	if (name == "alpha") {
	this->setAlpha(value);
	std::cout << "INFO in LauBelleSymNR::setResonanceParameter : Setting parameter alpha = " << this->getAlpha() << std::endl;
	} else {
	std::cerr << "WARNING in LauBelleSymNR::setResonanceParameter : Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauBelleSymNR::floatResonanceParameter(const TString& name)
	{
	if (name == "alpha") {
	if ( alpha_->fixed() ) {
	alpha_->fixed( kFALSE );
	this->addFloatingParameter( alpha_ );
	} else {
	std::cerr << "WARNING in LauBelleSymNR::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauBelleSymNR::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauBelleSymNR::getResonanceParameter(const TString& name)
	{
	if (name == "alpha") {
	return alpha_;
	} else {
	std::cerr << "WARNING in LauBelleSymNR::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauBelleSymNR::setAlpha(const Double_t alpha)
	{
	alpha_->value( alpha );
	alpha_->genValue( alpha );
	alpha_->initValue( alpha );
	}

	diff --git a/src/LauDabbaRes.cc b/src/LauDabbaRes.cc
	index a8f1341..743ce87 100644
	--- a/src/LauDabbaRes.cc
	+++ b/src/LauDabbaRes.cc
	@@ -1,250 +1,253 @@

	// Copyright University of Warwick 2010 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauDabbaRes.cc
	\brief File containing implementation of LauDabbaRes class.
	Formulae and data values from arXiv:0901.2217 - author D.V.Bugg
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauDabbaRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauDabbaRes)


	LauDabbaRes::LauDabbaRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	mSumSq_(0.0),
	sAdler_(0.0),
	b_(0),
	alpha_(0),
	beta_(0)
	{
	// Default constant factors
	const Double_t bVal = 24.49;
	const Double_t alphaVal = 0.1;
	const Double_t betaVal = 0.1;

	const TString& parNameBase = this->getSanitisedName();

	TString bName(parNameBase);
	bName += "_b";
	b_ = resInfo->getExtraParameter( bName );
	if ( b_ == 0 ) {
	b_ = new LauParameter( bName, bVal, 0.0, 100.0, kTRUE );
	+ b_->secondStage(kTRUE);
	resInfo->addExtraParameter( b_ );
	}

	TString alphaName(parNameBase);
	alphaName += "_alpha";
	alpha_ = resInfo->getExtraParameter( alphaName );
	if ( alpha_ == 0 ) {
	alpha_ = new LauParameter( alphaName, alphaVal, 0.0, 10.0, kTRUE );
	+ alpha_->secondStage(kTRUE);
	resInfo->addExtraParameter( alpha_ );
	}

	TString betaName(parNameBase);
	betaName += "_beta";
	beta_ = resInfo->getExtraParameter( betaName );
	if ( beta_ == 0 ) {
	beta_ = new LauParameter( betaName, betaVal, 0.0, 10.0, kTRUE );
	+ beta_->secondStage(kTRUE);
	resInfo->addExtraParameter( beta_ );
	}
	}

	LauDabbaRes::~LauDabbaRes()
	{
	}

	void LauDabbaRes::initialise()
	{
	// check that we have a D and a pi
	this->checkDaughterTypes();

	// Initialise various constants
	Double_t massDaug1 = this->getMassDaug1();
	Double_t massDaug2 = this->getMassDaug2();

	Double_t mSum = massDaug1 + massDaug2;
	mSumSq_ = mSum*mSum;

	Double_t massDaug1Sq = massDaug1*massDaug1;
	Double_t massDaug2Sq = massDaug2*massDaug2;
	sAdler_ = TMath::Max(massDaug1Sq,massDaug2Sq) - 0.5TMath::Min(massDaug1Sq,massDaug2Sq); // Adler zero at (mD)^2 - 0.5(mpi)^2

	Int_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauDabbaRes::initialise : Spin = " << resSpin << " is not zero! It will be ignored anyway!" << std::endl;
	}
	}

	void LauDabbaRes::checkDaughterTypes() const
	{
	// Check that the daughter tracks are D and pi. Otherwise issue a warning.
	Int_t resPairAmpInt = this->getPairInt();
	if (resPairAmpInt < 1 \|\| resPairAmpInt > 3) {
	std::cerr << "WARNING in LauDabbaRes::checkDaughterTypes : resPairAmpInt = " << resPairAmpInt << " is out of the range [1,2,3]." << std::endl;
	return;
	}

	// Check that daughter types agree
	const TString& nameDaug1 = this->getNameDaug1();
	const TString& nameDaug2 = this->getNameDaug2();
	if ( !( nameDaug1.Contains("pi", TString::kIgnoreCase) && nameDaug2.Contains("d", TString::kIgnoreCase) ) ) {
	if ( !( nameDaug2.Contains("pi", TString::kIgnoreCase) && nameDaug1.Contains("d", TString::kIgnoreCase) ) ) {
	std::cerr << "ERROR in LauDabbaRes::checkDaughterTypes : Dabba model is using daughters \"" << nameDaug1 << "\" and \"" << nameDaug2 << "\" that are not a D meson and a pion." << std::endl;
	}
	}
	}

	LauComplex LauDabbaRes::resAmp(Double_t mass, Double_t spinTerm)
	{
	// This function returns the complex dynamical amplitude for a Dabba distribution
	// given the invariant mass and cos(helicity) values.

	// Invariant mass squared combination for the system
	Double_t s = mass*mass;

	// Dabba is spin zero - so there are no helicity factors.
	// Just set it to 1.0 in case anyone decides to use it at a later date.
	spinTerm = 1.0;

	// Phase-space factor
	Double_t rho(0.0);
	Double_t sDiff = s - mSumSq_;
	if ( sDiff > 0.0 ) {
	rho = TMath::Sqrt(1.0 - mSumSq_/s);
	}

	const Double_t bVal = this->getBValue();
	const Double_t alphaVal = this->getAlphaValue();
	const Double_t betaVal = this->getBetaValue();

	Double_t realPart = 1.0 - betaVal * sDiff;

	Double_t imagPart = bVal * TMath::Exp( - alphaVal * sDiff ) * ( s - sAdler_ ) * rho;

	LauComplex resAmplitude( realPart, imagPart );

	Double_t denomFactor = realPartrealPart + imagPartimagPart;

	Double_t invDenomFactor = 0.0;
	if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}

	resAmplitude.rescale(spinTerm*invDenomFactor);

	return resAmplitude;
	}

	const std::vector<LauParameter*>& LauDabbaRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixBValue() ) {
	this->addFloatingParameter( b_ );
	}

	if ( ! this->fixAlphaValue() ) {
	this->addFloatingParameter( alpha_ );
	}

	if ( ! this->fixBetaValue() ) {
	this->addFloatingParameter( beta_ );
	}

	return this->getParameters();
	}

	void LauDabbaRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the lineshape
	if (name == "b") {
	this->setBValue(value);
	std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter b = " << this->getBValue() << std::endl;
	}
	else if (name == "alpha") {
	this->setAlphaValue(value);
	std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter alpha = " << this->getAlphaValue() << std::endl;
	}
	else if (name == "beta") {
	this->setBetaValue(value);
	std::cout << "INFO in LauDabbaRes::setResonanceParameter : Setting parameter beta = " << this->getBetaValue() << std::endl;
	}
	else {
	std::cerr << "WARNING in LauDabbaRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauDabbaRes::floatResonanceParameter(const TString& name)
	{
	if (name == "b") {
	if ( b_->fixed() ) {
	b_->fixed( kFALSE );
	this->addFloatingParameter( b_ );
	} else {
	std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "alpha") {
	if ( alpha_->fixed() ) {
	alpha_->fixed( kFALSE );
	this->addFloatingParameter( alpha_ );
	} else {
	std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "beta") {
	if ( beta_->fixed() ) {
	beta_->fixed( kFALSE );
	this->addFloatingParameter( beta_ );
	} else {
	std::cerr << "WARNING in LauDabbaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauDabbaRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauDabbaRes::getResonanceParameter(const TString& name)
	{
	if (name == "b") {
	return b_;
	} else if (name == "alpha") {
	return alpha_;
	} else if (name == "beta") {
	return beta_;
	} else {
	std::cerr << "WARNING in LauDabbaRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauDabbaRes::setBValue(const Double_t b)
	{
	b_->value( b );
	b_->genValue( b );
	b_->initValue( b );
	}

	void LauDabbaRes::setAlphaValue(const Double_t alpha)
	{
	alpha_->value( alpha );
	alpha_->genValue( alpha );
	alpha_->initValue( alpha );
	}

	void LauDabbaRes::setBetaValue(const Double_t beta)
	{
	beta_->value( beta );
	beta_->genValue( beta );
	beta_->initValue( beta );
	}

	diff --git a/src/LauFlatteRes.cc b/src/LauFlatteRes.cc
	index e86bccb..a123476 100644
	--- a/src/LauFlatteRes.cc
	+++ b/src/LauFlatteRes.cc
	@@ -1,281 +1,283 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauFlatteRes.cc
	\brief File containing implementation of LauFlatteRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauFlatteRes.hh"
	#include "LauResonanceInfo.hh"

	#include "TSystem.h"

	ClassImp(LauFlatteRes)


	LauFlatteRes::LauFlatteRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	g1_(0),
	g2_(0),
	mSumSq0_(0.0),
	mSumSq1_(0.0),
	mSumSq2_(0.0),
	mSumSq3_(0.0),
	useAdlerTerm_(kFALSE),
	sA_(0.0)
	{
	Double_t g1Val(0.);
	Double_t g2Val(0.);

	const TString& resName = this->getResonanceName();

	if ( resName == "f_0(980)" ) {
	mSumSq0_ = (LauConstants::mPi0 + LauConstants::mPi0) * (LauConstants::mPi0 + LauConstants::mPi0);
	mSumSq1_ = (LauConstants::mPi + LauConstants::mPi) * (LauConstants::mPi + LauConstants::mPi);
	mSumSq2_ = (LauConstants::mK + LauConstants::mK) * (LauConstants::mK + LauConstants::mK);
	mSumSq3_ = (LauConstants::mK0 + LauConstants::mK0) * (LauConstants::mK0 + LauConstants::mK0);
	// constant factors from BES data
	// resMass should be 0.965 +/- 0.008 +/- 0.006 GeV/c^2
	g1Val = 0.165; // +/- 0.010 +/- 0.015 GeV/c^2
	g2Val = g1Val*4.21; // +/- 0.25 +/- 0.21

	// or from E791
	//g1Val = 0.09;
	//g2Val = 0.02;

	// or from CERN/WA76
	//g1Val = 0.28;
	//g2Val = 0.56;
	} else if ( resName == "K*0_0(1430)" ) {
	mSumSq0_ = (LauConstants::mK0 + LauConstants::mPi0) * (LauConstants::mK0 + LauConstants::mPi0);
	mSumSq1_ = (LauConstants::mK + LauConstants::mPi) * (LauConstants::mK + LauConstants::mPi);
	mSumSq2_ = (LauConstants::mK0 + LauConstants::mEtaPrime) * (LauConstants::mK0 + LauConstants::mEtaPrime);
	mSumSq3_ = (LauConstants::mK0 + LauConstants::mEtaPrime) * (LauConstants::mK0 + LauConstants::mEtaPrime);
	//Phys. Lett. B 572, 1 (2003)
	// resMass should be 1.513 GeV/c^2
	g1Val = 0.304;
	g2Val = 0.380;
	useAdlerTerm_ = kTRUE;
	sA_ = 0.234;
	} else if ( resName == "K+_0(1430)" \|\| resName == "K-_0(1430)" ) {
	mSumSq0_ = (LauConstants::mK + LauConstants::mPi0) * (LauConstants::mK + LauConstants::mPi0);
	mSumSq1_ = (LauConstants::mK0 + LauConstants::mPi) * (LauConstants::mK0 + LauConstants::mPi);
	mSumSq2_ = (LauConstants::mK + LauConstants::mEtaPrime) * (LauConstants::mK + LauConstants::mEtaPrime);
	mSumSq3_ = (LauConstants::mK + LauConstants::mEtaPrime) * (LauConstants::mK + LauConstants::mEtaPrime);
	//Phys. Lett. B 572, 1 (2003)
	// resMass should be 1.513 GeV/c^2
	g1Val = 0.304;
	g2Val = 0.380;
	useAdlerTerm_ = kTRUE;
	sA_ = 0.234;
	} else if ( resName == "a0_0(980)" ) {
	mSumSq0_ = (LauConstants::mEta + LauConstants::mPi0) * (LauConstants::mEta + LauConstants::mPi0);
	mSumSq1_ = (LauConstants::mEta + LauConstants::mPi0) * (LauConstants::mEta + LauConstants::mPi0);
	mSumSq2_ = (LauConstants::mK + LauConstants::mK) * (LauConstants::mK + LauConstants::mK);
	mSumSq3_ = (LauConstants::mK0 + LauConstants::mK0) * (LauConstants::mK0 + LauConstants::mK0);
	//Phys. Rev. D 57, 3860 (1998)
	// resMass should be 0.982 +/- 0.003 GeV/c^2
	g1Val = 0.353;
	g2Val = g1Val*1.03;
	} else if ( resName == "a+_0(980)" \|\| resName == "a-_0(980)" ) {
	mSumSq0_ = (LauConstants::mEta + LauConstants::mPi) * (LauConstants::mEta + LauConstants::mPi);
	mSumSq1_ = (LauConstants::mEta + LauConstants::mPi) * (LauConstants::mEta + LauConstants::mPi);
	mSumSq2_ = (LauConstants::mK + LauConstants::mK0) * (LauConstants::mK + LauConstants::mK0);
	mSumSq3_ = (LauConstants::mK + LauConstants::mK0) * (LauConstants::mK + LauConstants::mK0);
	//Phys. Rev. D 57, 3860 (1998)
	g1Val = 0.353;
	g2Val = g1Val*1.03;
	}

	const TString& parNameBase = this->getSanitisedName();

	TString g1Name(parNameBase);
	g1Name += "_g1";
	g1_ = resInfo->getExtraParameter( g1Name );
	if ( g1_ == 0 ) {
	g1_ = new LauParameter( g1Name, g1Val, 0.0, 10.0, kTRUE );
	+ g1_->secondStage(kTRUE);
	resInfo->addExtraParameter( g1_ );
	}

	TString g2Name(parNameBase);
	g2Name += "_g2";
	g2_ = resInfo->getExtraParameter( g2Name );
	if ( g2_ == 0 ) {
	g2_ = new LauParameter( g2Name, g2Val, 0.0, 10.0, kTRUE );
	+ g2_->secondStage(kTRUE);
	resInfo->addExtraParameter( g2_ );
	}
	}

	LauFlatteRes::~LauFlatteRes()
	{
	delete g1_;
	delete g2_;
	}

	void LauFlatteRes::initialise()
	{
	const TString& resName = this->getResonanceName();
	if ( resName != "f_0(980)" && resName != "K0_0(1430)" && resName != "K+_0(1430)" && resName != "K*-_0(1430)" &&
	resName != "a0_0(980)" && resName != "a+_0(980)" && resName != "a-_0(980)" ) {
	std::cerr << "ERROR in LauFlatteRes::initialise : Unexpected resonance name \"" << resName << "\" for Flatte shape." << std::endl;
	gSystem->Exit(EXIT_FAILURE);
	}

	Int_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauFlatteRes::amplitude : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : Flatte amplitude is only defined for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}
	}

	LauComplex LauFlatteRes::resAmp(Double_t mass, Double_t spinTerm)
	{
	// This function returns the complex dynamical amplitude for a Flatte distribution
	// given the invariant mass and cos(helicity) values.

	const Double_t resMass = this->getMass();
	const Double_t resMassSq = resMass*resMass;
	const Double_t s = mass*mass; // Invariant mass squared combination for the system

	const Double_t g1Val = this->getg1Parameter();
	const Double_t g2Val = this->getg2Parameter();

	Double_t dMSq = resMassSq - s;
	Double_t rho1(0.0), rho2(0.0);
	if (s > mSumSq0_) {
	rho1 = TMath::Sqrt(1.0 - mSumSq0_/s)/3.0;
	if (s > mSumSq1_) {
	rho1 += 2.0*TMath::Sqrt(1.0 - mSumSq1_/s)/3.0;
	if (s > mSumSq2_) {
	rho2 = 0.5*TMath::Sqrt(1.0 - mSumSq2_/s);
	if (s > mSumSq3_) {
	rho2 += 0.5*TMath::Sqrt(1.0 - mSumSq3_/s);
	} else {
	// Continue analytically below higher channel thresholds
	// This contributes to the real part of the amplitude denominator
	dMSq += g2ValresMass0.5*TMath::Sqrt(mSumSq3_/s - 1.0);
	}
	} else {
	// Continue analytically below higher channel thresholds
	// This contributes to the real part of the amplitude denominator
	rho2 = 0.0;
	dMSq += g2ValresMass(0.5TMath::Sqrt(mSumSq2_/s - 1.0) + 0.5TMath::Sqrt(mSumSq3_/s - 1.0));
	}
	} else {
	// Continue analytically below higher channel thresholds
	// This contributes to the real part of the amplitude denominator
	dMSq += g1ValresMass2.0*TMath::Sqrt(mSumSq1_/s - 1.0)/3.0;
	}
	}

	Double_t massFactor = resMass;
	if(useAdlerTerm_) massFactor = ( s - sA_ ) / ( resMassSq - sA_ );
	const Double_t width1 = g1Valrho1massFactor;
	const Double_t width2 = g2Valrho2massFactor;
	const Double_t widthTerm = width1 + width2;

	LauComplex resAmplitude(dMSq, widthTerm);

	const Double_t denomFactor = dMSqdMSq + widthTermwidthTerm;

	Double_t invDenomFactor = 0.0;
	if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}

	resAmplitude.rescale(spinTerm*invDenomFactor);

	return resAmplitude;
	}

	const std::vector<LauParameter*>& LauFlatteRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixMass() ) {
	this->addFloatingParameter( this->getMassPar() );
	}

	// NB the width is given in terms of g1 and g2 so the normal width
	// parameter should be ignored, hence it is not added to the list

	if ( ! this->fixg1Parameter() ) {
	this->addFloatingParameter( g1_ );
	}

	if ( ! this->fixg2Parameter() ) {
	this->addFloatingParameter( g2_ );
	}

	return this->getParameters();
	}

	void LauFlatteRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	if (name == "g1") {
	this->setg1Parameter(value);
	std::cout << "INFO in LauFlatteRes::setResonanceParameter : Setting g1 parameter to " << this->getg1Parameter() << std::endl;
	} else if (name == "g2") {
	this->setg2Parameter(value);
	std::cout << "INFO in LauFlatteRes::setResonanceParameter : Setting g2 parameter to " << this->getg2Parameter() << std::endl;
	} else {
	std::cerr << "WARNING in LauFlatteRes::setResonanceParameter : Parameter name \"" << name << "\" not recognised." << std::endl;
	}
	}

	void LauFlatteRes::floatResonanceParameter(const TString& name)
	{
	if (name == "g1") {
	if ( g1_->fixed() ) {
	g1_->fixed( kFALSE );
	this->addFloatingParameter( g1_ );
	} else {
	std::cerr << "WARNING in LauFlatteRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "g2") {
	if ( g2_->fixed() ) {
	g2_->fixed( kFALSE );
	this->addFloatingParameter( g2_ );
	} else {
	std::cerr << "WARNING in LauFlatteRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauFlatteRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauFlatteRes::getResonanceParameter(const TString& name)
	{
	if (name == "g1") {
	return g1_;
	} else if (name == "g2") {
	return g2_;
	} else {
	std::cerr << "WARNING in LauFlatteRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauFlatteRes::setg1Parameter(const Double_t g1)
	{
	g1_->value( g1 );
	g1_->genValue( g1 );
	g1_->initValue( g1 );
	}

	void LauFlatteRes::setg2Parameter(const Double_t g2)
	{
	g2_->value( g2 );
	g2_->genValue( g2 );
	g2_->initValue( g2 );
	}

	diff --git a/src/LauKappaRes.cc b/src/LauKappaRes.cc
	index 0be066e..84702c3 100644
	--- a/src/LauKappaRes.cc
	+++ b/src/LauKappaRes.cc
	@@ -1,282 +1,286 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauKappaRes.cc
	\brief File containing implementation of LauKappaRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauKappaRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauKappaRes)


	LauKappaRes::LauKappaRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	mSumSq_((LauConstants::mPi+LauConstants::mK)*(LauConstants::mPi+LauConstants::mK)),
	sAdler_(LauConstants::mKSq-0.5*LauConstants::mPiSq),
	b1_(0),
	b2_(0),
	a_(0),
	m0_(0)
	{
	// Default constant factors from BES data
	const Double_t b1Val = 24.49;
	const Double_t b2Val = 0.0;
	const Double_t aVal = 2.5;
	const Double_t m0Val = 3.3;

	const TString& parNameBase = this->getSanitisedName();

	TString b1Name(parNameBase);
	b1Name += "_b1";
	b1_ = resInfo->getExtraParameter( b1Name );
	if ( b1_ == 0 ) {
	b1_ = new LauParameter( b1Name, b1Val, 0.0, 100.0, kTRUE );
	+ b1_->secondStage(kTRUE);
	resInfo->addExtraParameter( b1_ );
	}

	TString b2Name(parNameBase);
	b2Name += "_b2";
	b2_ = resInfo->getExtraParameter( b2Name );
	if ( b2_ == 0 ) {
	b2_ = new LauParameter( b2Name, b2Val, 0.0, 100.0, kTRUE );
	+ b2_->secondStage(kTRUE);
	resInfo->addExtraParameter( b2_ );
	}

	TString aName(parNameBase);
	aName += "_A";
	a_ = resInfo->getExtraParameter( aName );
	if ( a_ == 0 ) {
	a_ = new LauParameter( aName, aVal, 0.0, 10.0, kTRUE );
	+ a_->secondStage(kTRUE);
	resInfo->addExtraParameter( a_ );
	}

	TString m0Name(parNameBase);
	m0Name += "_m0";
	m0_ = resInfo->getExtraParameter( m0Name );
	if ( m0_ == 0 ) {
	m0_ = new LauParameter( m0Name, m0Val, 0.0, 10.0, kTRUE );
	+ m0_->secondStage(kTRUE);
	resInfo->addExtraParameter( m0_ );
	}
	}

	LauKappaRes::~LauKappaRes()
	{
	}

	void LauKappaRes::initialise()
	{
	this->checkDaughterTypes();

	Double_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauKappaRes::initialise : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : Kappa amplitude is only for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}
	}

	void LauKappaRes::checkDaughterTypes() const
	{
	// Check that the daughter tracks are K and pi. Otherwise issue a warning.

	Int_t resPairAmpInt = this->getPairInt();
	if (resPairAmpInt < 1 \|\| resPairAmpInt > 3) {
	std::cerr << "WARNING in LauKappaRes::checkDaughterTypes : resPairAmpInt = " << resPairAmpInt << " is out of the range [1,2,3]." << std::endl;
	return;
	}

	const TString& nameDaug1 = this->getNameDaug1();
	const TString& nameDaug2 = this->getNameDaug2();
	if ( !( nameDaug1.Contains("pi", TString::kIgnoreCase) && nameDaug2.Contains("k", TString::kIgnoreCase) ) ) {
	if ( !( nameDaug2.Contains("pi", TString::kIgnoreCase) && nameDaug1.Contains("k", TString::kIgnoreCase) ) ) {
	std::cerr << "ERROR in LauKappaRes::checkDaughterTypes : Kappa model is using daughters \"" << nameDaug1 << "\" and \"" << nameDaug2 << "\" that are not a kaon and a pion." << std::endl;
	}
	}
	}

	LauComplex LauKappaRes::resAmp(Double_t mass, Double_t spinTerm)
	{
	// This function returns the complex dynamical amplitude for a Kappa distribution
	// given the invariant mass and cos(helicity) values.

	// First check that the appropriate daughters are either pi+pi- or K+K-
	// Check that the daughter tracks are the same type. Otherwise issue a warning
	// and set the type to be pion for the Kappa distribution. Returns the
	// integer defined by the enum LauKappaRes::KappaPartType.

	Double_t s = mass*mass; // Invariant mass squared combination for the system
	Double_t rho(0.0); // Phase-space factor
	if (s > mSumSq_) {rho = TMath::Sqrt(1.0 - mSumSq_/s);}

	const Double_t m0Val = this->getM0Value();
	const Double_t m0Sq = m0Val * m0Val;

	const Double_t aVal = this->getAValue();
	const Double_t b1Val = this->getB1Value();
	const Double_t b2Val = this->getB2Value();

	Double_t f = b2Val*s + b1Val; // f(s) function
	Double_t numerator = s - sAdler_;
	Double_t denom = m0Sq - sAdler_;
	Double_t gamma(0.0);
	if (TMath::Abs(denom) > 1e-10 && TMath::Abs(aVal) > 1e-10) {
	// Decay width of the system
	gamma = rho(numerator/denom)f*TMath::Exp(-(s - m0Sq)/aVal);
	}

	// Now form the complex amplitude - use relativistic BW form (without barrier factors)
	// Note that the M factor in the denominator is not the "pole" at ~500 MeV, but is
	// m0_ = 0.9264, the mass when the phase shift goes through 90 degrees.

	Double_t dMSq = m0Sq - s;
	Double_t widthTerm = gamma*m0Val;
	LauComplex resAmplitude(dMSq, widthTerm);

	Double_t denomFactor = dMSqdMSq + widthTermwidthTerm;

	Double_t invDenomFactor = 0.0;
	if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}

	resAmplitude.rescale(spinTerm*invDenomFactor);

	return resAmplitude;
	}

	const std::vector<LauParameter*>& LauKappaRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixB1Value() ) {
	this->addFloatingParameter( b1_ );
	}

	if ( ! this->fixB2Value() ) {
	this->addFloatingParameter( b2_ );
	}

	if ( ! this->fixAValue() ) {
	this->addFloatingParameter( a_ );
	}

	if ( ! this->fixM0Value() ) {
	this->addFloatingParameter( m0_ );
	}

	return this->getParameters();
	}

	void LauKappaRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the lineshape
	if (name == "b1") {
	this->setB1Value(value);
	std::cout << "INFO in LauKappaRes::setResonanceParameter : Setting parameter b1 = " << this->getB1Value() << std::endl;
	}
	else if (name == "b2") {
	this->setB2Value(value);
	std::cout << "INFO in LauKappaRes::setResonanceParameter : Setting parameter b2 = " << this->getB2Value() << std::endl;
	}
	else if (name == "A") {
	this->setAValue(value);
	std::cout << "INFO in LauKappaRes::setResonanceParameter : Setting parameter A = " << this->getAValue() << std::endl;
	}
	else if (name == "m0") {
	this->setM0Value(value);
	std::cout << "INFO in LauKappaRes::setResonanceParameter : Setting parameter m0 = " << this->getM0Value() << std::endl;
	}
	else {
	std::cerr << "WARNING in LauKappaRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauKappaRes::floatResonanceParameter(const TString& name)
	{
	if (name == "b1") {
	if ( b1_->fixed() ) {
	b1_->fixed( kFALSE );
	this->addFloatingParameter( b1_ );
	} else {
	std::cerr << "WARNING in LauKappaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "b2") {
	if ( b2_->fixed() ) {
	b2_->fixed( kFALSE );
	this->addFloatingParameter( b2_ );
	} else {
	std::cerr << "WARNING in LauKappaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "A") {
	if ( a_->fixed() ) {
	a_->fixed( kFALSE );
	this->addFloatingParameter( a_ );
	} else {
	std::cerr << "WARNING in LauKappaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "m0") {
	if ( m0_->fixed() ) {
	m0_->fixed( kFALSE );
	this->addFloatingParameter( m0_ );
	} else {
	std::cerr << "WARNING in LauKappaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauKappaRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauKappaRes::getResonanceParameter(const TString& name)
	{
	if (name == "b1") {
	return b1_;
	} else if (name == "b2") {
	return b2_;
	} else if (name == "A") {
	return a_;
	} else if (name == "m0") {
	return m0_;
	} else {
	std::cerr << "WARNING in LauKappaRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauKappaRes::setB1Value(const Double_t b1)
	{
	b1_->value( b1 );
	b1_->genValue( b1 );
	b1_->initValue( b1 );
	}

	void LauKappaRes::setB2Value(const Double_t b2)
	{
	b2_->value( b2 );
	b2_->genValue( b2 );
	b2_->initValue( b2 );
	}

	void LauKappaRes::setAValue(const Double_t A)
	{
	a_->value( A );
	a_->genValue( A );
	a_->initValue( A );
	}

	void LauKappaRes::setM0Value(const Double_t m0)
	{
	m0_->value( m0 );
	m0_->genValue( m0 );
	m0_->initValue( m0 );
	}

	diff --git a/src/LauLASSBWRes.cc b/src/LauLASSBWRes.cc
	index 8d05a5e..d295ada 100644
	--- a/src/LauLASSBWRes.cc
	+++ b/src/LauLASSBWRes.cc
	@@ -1,228 +1,230 @@

	// Copyright University of Warwick 2008 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauLASSBWRes.cc
	\brief File containing implementation of LauLASSBWRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauLASSBWRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauLASSBWRes)


	LauLASSBWRes::LauLASSBWRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	q0_(0.0),
	mDaugSum_(0.0),
	mDaugSumSq_(0.0),
	mDaugDiff_(0.0),
	mDaugDiffSq_(0.0),
	resMassSq_(0.0),
	r_(0),
	a_(0)
	{
	// Default values for LASS parameters
	const Double_t rVal = 3.32;
	const Double_t aVal = 2.07;

	const TString& parNameBase = this->getSanitisedName();

	TString rName(parNameBase);
	rName += "_r";
	r_ = resInfo->getExtraParameter( rName );
	if ( r_ == 0 ) {
	r_ = new LauParameter( rName, rVal, 0.0, 10.0, kTRUE );
	+ r_->secondStage(kTRUE);
	resInfo->addExtraParameter( r_ );
	}

	TString aName(parNameBase);
	aName += "_a";
	a_ = resInfo->getExtraParameter( aName );
	if ( a_ == 0 ) {
	a_ = new LauParameter( aName, aVal, 0.0, 10.0, kTRUE );
	+ a_->secondStage(kTRUE);
	resInfo->addExtraParameter( a_ );
	}
	}

	LauLASSBWRes::~LauLASSBWRes()
	{
	delete r_;
	delete a_;
	}

	void LauLASSBWRes::initialise()
	{
	// Create the mass sums and differences
	Double_t massDaug1 = this->getMassDaug1();
	Double_t massDaug2 = this->getMassDaug2();

	mDaugSum_ = massDaug1 + massDaug2;
	mDaugSumSq_ = mDaugSum_*mDaugSum_;

	mDaugDiff_ = massDaug1 - massDaug2;
	mDaugDiffSq_ = mDaugDiff_*mDaugDiff_;

	Int_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauLASSBWRes::amplitude : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : LASS amplitude is only for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}

	this->calcQ0();
	}

	void LauLASSBWRes::calcQ0()
	{
	// Decay momentum of either daughter in the resonance rest frame
	// when resonance mass = rest-mass value, m_0 (PDG value)

	resMass_ = this->getMass();
	resMassSq_ = resMass_*resMass_;

	q0_ = TMath::Sqrt((resMassSq_ - mDaugSumSq_)(resMassSq_ - mDaugDiffSq_))/(2.0resMass_);
	}

	LauComplex LauLASSBWRes::resAmp(Double_t mass, Double_t /spinTerm/)
	{
	LauComplex resAmplitude(0.0, 0.0);

	if (mass < 1e-10) {
	std::cerr << "WARNING in LauLASSBWRes::amplitude : mass < 1e-10." << std::endl;
	return LauComplex(0.0, 0.0);
	}

	// Calculate the width of the resonance (as a function of mass)
	// q is the momentum of either daughter in the resonance rest-frame
	const Double_t q = this->getQ();
	const Double_t resMass = this->getMass();
	const Double_t resWidth = this->getWidth();

	// If the mass is floating and their value have changed
	// we need to recalculate everything that assumes this value
	if ( (!this->fixMass()) && resMass != resMass_ ) {
	this->calcQ0();
	}

	Double_t qRatio = q/q0_;
	Double_t totWidth = resWidthqRatio(resMass/mass);

	Double_t massSq = mass*mass;
	Double_t massSqTerm = resMassSq_ - massSq;

	// Compute the complex amplitude
	resAmplitude = LauComplex(massSqTerm, resMass*totWidth);

	// Scale by the denominator factor
	resAmplitude.rescale((resMassSq_resWidth/q0_)/(massSqTermmassSqTerm + resMassSq_totWidthtotWidth));

	// Calculate the phase shift term
	const Double_t rVal = this->getEffectiveRange();
	const Double_t aVal = this->getScatteringLength();
	const Double_t tandeltaB = (2.0aValq)/(2.0 + aValrValq*q);
	const Double_t tanSq = tandeltaB*tandeltaB;
	const Double_t cos2PhaseShift = (1.0 - tanSq) / (1.0 + tanSq);
	const Double_t sin2PhaseShift = 2.0*tandeltaB / (1.0 + tanSq);
	LauComplex phaseShift(cos2PhaseShift, sin2PhaseShift);

	// Multiply by the phase shift
	resAmplitude = resAmplitude * phaseShift;

	return resAmplitude;
	}

	const std::vector<LauParameter*>& LauLASSBWRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixMass() ) {
	this->addFloatingParameter( this->getMassPar() );
	}

	if ( ! this->fixWidth() ) {
	this->addFloatingParameter( this->getWidthPar() );
	}

	if ( ! this->fixEffectiveRange() ) {
	this->addFloatingParameter( r_ );
	}

	if ( ! this->fixScatteringLength() ) {
	this->addFloatingParameter( a_ );
	}

	return this->getParameters();
	}

	void LauLASSBWRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the LASS lineshape dynamics
	if (name == "a") {
	this->setScatteringLength(value);
	std::cout << "INFO in LauLASSBWRes::setResonanceParameter : Setting LASS Scattering Length = " << this->getScatteringLength() << std::endl;
	} else if (name == "r") {
	this->setEffectiveRange(value);
	std::cout << "INFO in LauLASSBWRes::setResonanceParameter : Setting LASS Effective Range = " << this->getEffectiveRange() << std::endl;
	} else {
	std::cerr << "WARNING in LauLASSBWRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauLASSBWRes::floatResonanceParameter(const TString& name)
	{
	if (name == "a") {
	if ( a_->fixed() ) {
	a_->fixed( kFALSE );
	this->addFloatingParameter( a_ );
	} else {
	std::cerr << "WARNING in LauLASSBWRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "r") {
	if ( r_->fixed() ) {
	r_->fixed( kFALSE );
	this->addFloatingParameter( r_ );
	} else {
	std::cerr << "WARNING in LauLASSBWRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauLASSBWRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauLASSBWRes::getResonanceParameter(const TString& name)
	{
	if (name == "a") {
	return a_;
	} else if (name == "r") {
	return r_;
	} else {
	std::cerr << "WARNING in LauLASSBWRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauLASSBWRes::setEffectiveRange(const Double_t r)
	{
	r_->value( r );
	r_->genValue( r );
	r_->initValue( r );
	}

	void LauLASSBWRes::setScatteringLength(const Double_t a)
	{
	a_->value( a );
	a_->genValue( a );
	a_->initValue( a );
	}

	diff --git a/src/LauLASSNRRes.cc b/src/LauLASSNRRes.cc
	index 3795003..a49b593 100644
	--- a/src/LauLASSNRRes.cc
	+++ b/src/LauLASSNRRes.cc
	@@ -1,182 +1,184 @@

	// Copyright University of Warwick 2008 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauLASSNRRes.cc
	\brief File containing implementation of LauLASSNRRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauLASSNRRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauLASSNRRes)


	LauLASSNRRes::LauLASSNRRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	r_(0),
	a_(0),
	cutOff_(0.0)
	{
	// Default values for LASS parameters
	cutOff_ = 1.8;
	const Double_t rVal = 3.32;
	const Double_t aVal = 2.07;

	const TString& parNameBase = this->getSanitisedName();

	TString rName(parNameBase);
	rName += "_r";
	r_ = resInfo->getExtraParameter( rName );
	if ( r_ == 0 ) {
	r_ = new LauParameter( rName, rVal, 0.0, 10.0, kTRUE );
	+ r_->secondStage(kTRUE);
	resInfo->addExtraParameter( r_ );
	}

	TString aName(parNameBase);
	aName += "_a";
	a_ = resInfo->getExtraParameter( aName );
	if ( a_ == 0 ) {
	a_ = new LauParameter( aName, aVal, 0.0, 10.0, kTRUE );
	+ a_->secondStage(kTRUE);
	resInfo->addExtraParameter( a_ );
	}
	}

	LauLASSNRRes::~LauLASSNRRes()
	{
	delete r_;
	delete a_;
	}

	void LauLASSNRRes::initialise()
	{
	Int_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauLASSNRRes::amplitude : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : LASS amplitude is only for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}
	}

	LauComplex LauLASSNRRes::resAmp(Double_t mass, Double_t /spinTerm/)
	{
	LauComplex bkgAmplitude(0.0, 0.0);

	if (mass < 1e-10) {
	std::cerr << "WARNING in LauLASSNRRes::amplitude : mass < 1e-10." << std::endl;
	return LauComplex(0.0, 0.0);
	}

	if (mass < cutOff_) {
	// q is the momentum of either daughter in the resonance rest-frame
	const Double_t q = this->getQ();

	// Calculate the phase shift term
	const Double_t rVal = this->getEffectiveRange();
	const Double_t aVal = this->getScatteringLength();
	const Double_t qcotdeltaB = 1.0/aVal + (rValqq)/2.0;

	// Compute the complex amplitude
	bkgAmplitude = LauComplex(qcotdeltaB, q);

	// Scale by the numerator and denominator factors
	bkgAmplitude.rescale(mass/(qcotdeltaBqcotdeltaB + qq));
	}

	return bkgAmplitude;

	}

	const std::vector<LauParameter*>& LauLASSNRRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixMass() ) {
	this->addFloatingParameter( this->getMassPar() );
	}

	if ( ! this->fixWidth() ) {
	this->addFloatingParameter( this->getWidthPar() );
	}

	if ( ! this->fixEffectiveRange() ) {
	this->addFloatingParameter( r_ );
	}

	if ( ! this->fixScatteringLength() ) {
	this->addFloatingParameter( a_ );
	}

	return this->getParameters();
	}

	void LauLASSNRRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the LASS lineshape dynamics
	if (name == "a") {
	this->setScatteringLength(value);
	std::cout << "INFO in LauLASSNRRes::setResonanceParameter : Setting LASS Scattering Length = " << this->getScatteringLength() << std::endl;
	} else if (name == "r") {
	this->setEffectiveRange(value);
	std::cout << "INFO in LauLASSNRRes::setResonanceParameter : Setting LASS Effective Range = " << this->getEffectiveRange() << std::endl;
	} else {
	std::cerr << "WARNING in LauLASSNRRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauLASSNRRes::floatResonanceParameter(const TString& name)
	{
	if (name == "a") {
	if ( a_->fixed() ) {
	a_->fixed( kFALSE );
	this->addFloatingParameter( a_ );
	} else {
	std::cerr << "WARNING in LauLASSNRRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "r") {
	if ( r_->fixed() ) {
	r_->fixed( kFALSE );
	this->addFloatingParameter( r_ );
	} else {
	std::cerr << "WARNING in LauLASSNRRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauLASSNRRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauLASSNRRes::getResonanceParameter(const TString& name)
	{
	if (name == "a") {
	return a_;
	} else if (name == "r") {
	return r_;
	} else {
	std::cerr << "WARNING in LauLASSNRRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauLASSNRRes::setEffectiveRange(const Double_t r)
	{
	r_->value( r );
	r_->genValue( r );
	r_->initValue( r );
	}

	void LauLASSNRRes::setScatteringLength(const Double_t a)
	{
	a_->value( a );
	a_->genValue( a );
	a_->initValue( a );
	}

	diff --git a/src/LauLASSRes.cc b/src/LauLASSRes.cc
	index 6fc770f..dd794b7 100644
	--- a/src/LauLASSRes.cc
	+++ b/src/LauLASSRes.cc
	@@ -1,261 +1,263 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauLASSRes.cc
	\brief File containing implementation of LauLASSRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauLASSRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauLASSRes)


	LauLASSRes::LauLASSRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	q0_(0.0),
	mDaugSum_(0.0),
	mDaugSumSq_(0.0),
	mDaugDiff_(0.0),
	mDaugDiffSq_(0.0),
	resMassSq_(0.0),
	r_(0),
	a_(0),
	cutOff_(0.0)
	{
	// Default values for LASS parameters
	cutOff_ = 1.8;
	const Double_t rVal = 3.32;
	const Double_t aVal = 2.07;

	const TString& parNameBase = this->getSanitisedName();

	TString rName(parNameBase);
	rName += "_r";
	r_ = resInfo->getExtraParameter( rName );
	if ( r_ == 0 ) {
	r_ = new LauParameter( rName, rVal, 0.0, 10.0, kTRUE );
	+ r_->secondStage(kTRUE);
	resInfo->addExtraParameter( r_ );
	}

	TString aName(parNameBase);
	aName += "_a";
	a_ = resInfo->getExtraParameter( aName );
	if ( a_ == 0 ) {
	a_ = new LauParameter( aName, aVal, 0.0, 10.0, kTRUE );
	+ a_->secondStage(kTRUE);
	resInfo->addExtraParameter( a_ );
	}
	}

	LauLASSRes::~LauLASSRes()
	{
	delete r_;
	delete a_;
	}

	void LauLASSRes::initialise()
	{
	// Create the mass sums and differences
	Double_t massDaug1 = this->getMassDaug1();
	Double_t massDaug2 = this->getMassDaug2();

	mDaugSum_ = massDaug1 + massDaug2;
	mDaugSumSq_ = mDaugSum_*mDaugSum_;

	mDaugDiff_ = massDaug1 - massDaug2;
	mDaugDiffSq_ = mDaugDiff_*mDaugDiff_;

	Int_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauLASSRes::amplitude : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : LASS amplitude is only for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}

	this->calcQ0();
	}

	void LauLASSRes::calcQ0()
	{
	// Decay momentum of either daughter in the resonance rest frame
	// when resonance mass = rest-mass value, m_0 (PDG value)

	resMass_ = this->getMass();
	resMassSq_ = resMass_*resMass_;

	q0_ = TMath::Sqrt((resMassSq_ - mDaugSumSq_)(resMassSq_ - mDaugDiffSq_))/(2.0resMass_);
	}

	LauComplex LauLASSRes::resAmp(Double_t mass, Double_t /spinTerm/)
	{
	LauComplex resAmplitude(0.0, 0.0);
	LauComplex bkgAmplitude(0.0, 0.0);
	LauComplex totAmplitude(0.0, 0.0);

	if (mass < 1e-10) {
	std::cerr << "WARNING in LauLASSRes::amplitude : mass < 1e-10." << std::endl;
	return LauComplex(0.0, 0.0);
	}

	//---------------------------
	// First do the resonant part
	//---------------------------

	// Calculate the width of the resonance (as a function of mass)
	// q is the momentum of either daughter in the resonance rest-frame
	const Double_t q = this->getQ();
	const Double_t resMass = this->getMass();
	const Double_t resWidth = this->getWidth();

	// If the mass is floating and their value have changed
	// we need to recalculate everything that assumes this value
	if ( (!this->fixMass()) && resMass != resMass_ ) {
	this->calcQ0();
	}

	Double_t qRatio = q/q0_;
	Double_t totWidth = resWidthqRatio(resMass/mass);

	Double_t massSq = mass*mass;
	Double_t massSqTerm = resMassSq_ - massSq;

	// Compute the complex amplitude
	resAmplitude = LauComplex(massSqTerm, resMass*totWidth);

	// Scale by the denominator factor
	resAmplitude.rescale((resMassSq_resWidth/q0_)/(massSqTermmassSqTerm + resMassSq_totWidthtotWidth));

	// Calculate the phase shift term
	const Double_t rVal = this->getEffectiveRange();
	const Double_t aVal = this->getScatteringLength();
	const Double_t tandeltaB = (2.0aValq)/(2.0 + aValrValq*q);
	const Double_t tanSq = tandeltaB*tandeltaB;
	const Double_t cos2PhaseShift = (1.0 - tanSq) / (1.0 + tanSq);
	const Double_t sin2PhaseShift = 2.0*tandeltaB / (1.0 + tanSq);
	LauComplex phaseShift(cos2PhaseShift, sin2PhaseShift);

	// Multiply by the phase shift
	resAmplitude = resAmplitude * phaseShift;


	//--------------------------------
	// Now do the effective range part
	//--------------------------------

	// Form the real and imaginary parts
	const Double_t qcotdeltaB = 1.0/aVal + (rValqq)/2.0;

	// Compute the complex amplitude
	bkgAmplitude = LauComplex(qcotdeltaB, q);

	// Scale by the numerator and denominator factors
	bkgAmplitude.rescale(mass/(qcotdeltaBqcotdeltaB + qq));


	//------------------
	// Add them together
	//------------------

	if (mass > cutOff_) {
	totAmplitude = resAmplitude;
	} else {
	totAmplitude = bkgAmplitude + resAmplitude;
	}

	return totAmplitude;
	}

	const std::vector<LauParameter*>& LauLASSRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixMass() ) {
	this->addFloatingParameter( this->getMassPar() );
	}

	if ( ! this->fixWidth() ) {
	this->addFloatingParameter( this->getWidthPar() );
	}

	if ( ! this->fixEffectiveRange() ) {
	this->addFloatingParameter( r_ );
	}

	if ( ! this->fixScatteringLength() ) {
	this->addFloatingParameter( a_ );
	}

	return this->getParameters();
	}

	void LauLASSRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the LASS lineshape dynamics
	if (name == "a") {
	this->setScatteringLength(value);
	std::cout << "INFO in LauLASSRes::setResonanceParameter : Setting LASS Scattering Length = " << this->getScatteringLength() << std::endl;
	} else if (name == "r") {
	this->setEffectiveRange(value);
	std::cout << "INFO in LauLASSRes::setResonanceParameter : Setting LASS Effective Range = " << this->getEffectiveRange() << std::endl;
	} else {
	std::cerr << "WARNING in LauLASSRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauLASSRes::floatResonanceParameter(const TString& name)
	{
	if (name == "a") {
	if ( a_->fixed() ) {
	a_->fixed( kFALSE );
	this->addFloatingParameter( a_ );
	} else {
	std::cerr << "WARNING in LauLASSRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "r") {
	if ( r_->fixed() ) {
	r_->fixed( kFALSE );
	this->addFloatingParameter( r_ );
	} else {
	std::cerr << "WARNING in LauLASSRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauLASSRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauLASSRes::getResonanceParameter(const TString& name)
	{
	if (name == "a") {
	return a_;
	} else if (name == "r") {
	return r_;
	} else {
	std::cerr << "WARNING in LauLASSRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauLASSRes::setEffectiveRange(const Double_t r)
	{
	r_->value( r );
	r_->genValue( r );
	r_->initValue( r );
	}

	void LauLASSRes::setScatteringLength(const Double_t a)
	{
	a_->value( a );
	a_->genValue( a );
	a_->initValue( a );
	}

	diff --git a/src/LauNRAmplitude.cc b/src/LauNRAmplitude.cc
	index 4da120d..e13f845 100644
	--- a/src/LauNRAmplitude.cc
	+++ b/src/LauNRAmplitude.cc
	@@ -1,271 +1,276 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauNRAmplitude.cc
	\brief File containing implementation of LauNRAmplitude class.
	*/

	#include <iostream>

	#include "LauKinematics.hh"
	#include "LauNRAmplitude.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauNRAmplitude)

	LauNRAmplitude::LauNRAmplitude(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	d_(0),
	c1_(0),
	c2_(0),
	p1_(0),
	p2_(0)
	{
	// Default values for parameters, taken from arXiv:0709.0075v1 [hep-ph]
	const Double_t dVal(1.3232e-3);
	const Double_t c1Val(0.65);
	const Double_t c2Val(0.55);
	const Double_t p1Val(18.0);
	const Double_t p2Val(15.0);

	const TString& parNameBase = this->getSanitisedName();

	TString dName(parNameBase);
	dName += "_d";
	d_ = resInfo->getExtraParameter( dName );
	if ( d_ == 0 ) {
	d_ = new LauParameter( dName, dVal, 0.0, 1.0, kTRUE );
	+ d_->secondStage(kTRUE);
	resInfo->addExtraParameter( d_ );
	}

	TString c1Name(parNameBase);
	c1Name += "_c1";
	c1_ = resInfo->getExtraParameter( c1Name );
	if ( c1_ == 0 ) {
	c1_ = new LauParameter( c1Name, c1Val, 0.0, 2.0, kTRUE );
	+ c1_->secondStage(kTRUE);
	resInfo->addExtraParameter( c1_ );
	}

	TString c2Name(parNameBase);
	c2Name += "_c2";
	c2_ = resInfo->getExtraParameter( c2Name );
	if ( c2_ == 0 ) {
	c2_ = new LauParameter( c2Name, c2Val, 0.0, 2.0, kTRUE );
	+ c2_->secondStage(kTRUE);
	resInfo->addExtraParameter( c2_ );
	}

	TString p1Name(parNameBase);
	p1Name += "_p1";
	p1_ = resInfo->getExtraParameter( p1Name );
	if ( p1_ == 0 ) {
	p1_ = new LauParameter( p1Name, p1Val, 0.0, 50.0, kTRUE );
	+ p1_->secondStage(kTRUE);
	resInfo->addExtraParameter( p1_ );
	}

	TString p2Name(parNameBase);
	p2Name += "_p2";
	p2_ = resInfo->getExtraParameter( p2Name );
	if ( p2_ == 0 ) {
	p2_ = new LauParameter( p2Name, p2Val, 0.0, 50.0, kTRUE );
	+ p2_->secondStage(kTRUE);
	resInfo->addExtraParameter( p2_ );
	}
	}

	LauNRAmplitude::~LauNRAmplitude()
	{
	delete d_;
	delete c1_;
	delete c2_;
	delete p1_;
	delete p2_;
	}

	void LauNRAmplitude::initialise()
	{
	}

	LauComplex LauNRAmplitude::amplitude(const LauKinematics* kinematics)
	{
	// Get the information from the kinematics object
	const Double_t m13Sq = kinematics->getm13Sq();
	const Double_t m23Sq = kinematics->getm23Sq();
	const Double_t m13 = kinematics->getm13();
	const Double_t m23 = kinematics->getm23();

	// Calculate the magnitude
	Double_t magnitude = TMath::Sqrt( m13 * m23 *
	this->f(m23Sq, c1_->value(), p1_->value()) *
	this->f(m13Sq, c2_->value(), p2_->value()) *
	TMath::Exp( -1.0 * d_->value() * m13Sqm13Sq m23Sq*m23Sq )
	);

	// return the amplitude
	LauComplex resAmplitude(magnitude, 0.0);
	return resAmplitude;
	}

	LauComplex LauNRAmplitude::resAmp(Double_t mass, Double_t spinTerm)
	{
	std::cerr << "ERROR in LauNRAmplitude::resAmp : This method shouldn't get called." << std::endl;
	std::cerr << " Returning zero amplitude for mass = " << mass << " and spinTerm = " << spinTerm << "." << std::endl;
	return LauComplex(0.0, 0.0);
	}

	Double_t LauNRAmplitude::f(const Double_t s, const Double_t c, const Double_t p) const
	{
	return 1.0 / (1.0 + TMath::Exp( c * (s-p) ));
	}

	// TODO up to here!
	const std::vector<LauParameter*>& LauNRAmplitude::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixdParameter() ) {
	this->addFloatingParameter( d_ );
	}

	if ( ! this->fixc1Parameter() ) {
	this->addFloatingParameter( c1_ );
	}

	if ( ! this->fixc2Parameter() ) {
	this->addFloatingParameter( c2_ );
	}

	if ( ! this->fixp1Parameter() ) {
	this->addFloatingParameter( p1_ );
	}

	if ( ! this->fixp2Parameter() ) {
	this->addFloatingParameter( p2_ );
	}

	return this->getParameters();
	}

	void LauNRAmplitude::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the NRAmplitude lineshape dynamics
	if (name == "d") {
	this->setdParameter(value);
	std::cout << "INFO in LauNRAmplitude::setResonanceParameter : Setting NRAmplitude d = " << this->getdParameter() << std::endl;
	} else if (name == "c1") {
	this->setc1Parameter(value);
	std::cout << "INFO in LauNRAmplitude::setResonanceParameter : Setting NRAmplitude c1 = " << this->getc1Parameter() << std::endl;
	} else if (name == "c2") {
	this->setc2Parameter(value);
	std::cout << "INFO in LauNRAmplitude::setResonanceParameter : Setting NRAmplitude c2 = " << this->getc2Parameter() << std::endl;
	} else if (name == "p1") {
	this->setp1Parameter(value);
	std::cout << "INFO in LauNRAmplitude::setResonanceParameter : Setting NRAmplitude p1 = " << this->getp1Parameter() << std::endl;
	} else if (name == "p2") {
	this->setp2Parameter(value);
	std::cout << "INFO in LauNRAmplitude::setResonanceParameter : Setting NRAmplitude p2 = " << this->getp2Parameter() << std::endl;
	} else {
	std::cerr << "WARNING in LauNRAmplitude::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauNRAmplitude::floatResonanceParameter(const TString& name)
	{
	if (name == "d") {
	if ( d_->fixed() ) {
	d_->fixed( kFALSE );
	this->addFloatingParameter( d_ );
	} else {
	std::cerr << "WARNING in LauNRAmplitude::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "c1") {
	if ( c1_->fixed() ) {
	c1_->fixed( kFALSE );
	this->addFloatingParameter( c1_ );
	} else {
	std::cerr << "WARNING in LauNRAmplitude::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "c2") {
	if ( c2_->fixed() ) {
	c2_->fixed( kFALSE );
	this->addFloatingParameter( c2_ );
	} else {
	std::cerr << "WARNING in LauNRAmplitude::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "p1") {
	if ( p1_->fixed() ) {
	p1_->fixed( kFALSE );
	this->addFloatingParameter( p1_ );
	} else {
	std::cerr << "WARNING in LauNRAmplitude::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "p2") {
	if ( p2_->fixed() ) {
	p2_->fixed( kFALSE );
	this->addFloatingParameter( p2_ );
	} else {
	std::cerr << "WARNING in LauNRAmplitude::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauNRAmplitude::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauNRAmplitude::getResonanceParameter(const TString& name)
	{
	if (name == "d") {
	return d_;
	} else if (name == "c1") {
	return c1_;
	} else if (name == "c2") {
	return c2_;
	} else if (name == "p1") {
	return p1_;
	} else if (name == "p2") {
	return p2_;
	} else {
	std::cerr << "WARNING in LauNRAmplitude::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauNRAmplitude::setdParameter(const Double_t d)
	{
	d_->value( d );
	d_->genValue( d );
	d_->initValue( d );
	}

	void LauNRAmplitude::setc1Parameter(const Double_t c1)
	{
	c1_->value( c1 );
	c1_->genValue( c1 );
	c1_->initValue( c1 );
	}

	void LauNRAmplitude::setc2Parameter(const Double_t c2)
	{
	c2_->value( c2 );
	c2_->genValue( c2 );
	c2_->initValue( c2 );
	}

	void LauNRAmplitude::setp1Parameter(const Double_t p1)
	{
	p1_->value( p1 );
	p1_->genValue( p1 );
	p1_->initValue( p1 );
	}

	void LauNRAmplitude::setp2Parameter(const Double_t p2)
	{
	p2_->value( p2 );
	p2_->genValue( p2 );
	p2_->initValue( p2 );
	}

	diff --git a/src/LauResonanceInfo.cc b/src/LauResonanceInfo.cc
	index ef26215..1bdc183 100644
	--- a/src/LauResonanceInfo.cc
	+++ b/src/LauResonanceInfo.cc
	@@ -1,184 +1,186 @@

	// Copyright University of Warwick 2006 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauResonanceInfo.cc
	\brief File containing implementation of LauResonanceInfo class.
	*/

	#include <iostream>

	#include "LauParameter.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauResonanceInfo)


	LauResonanceInfo::LauResonanceInfo(const TString& name, Double_t mass, Double_t width, Int_t spin, Int_t charge, Double_t range) :
	name_(name),
	sanitisedName_(""),
	mass_(0),
	width_(0),
	spin_(spin),
	charge_(charge),
	range_(range),
	conjugate_(0),
	extraPars_()
	{
	this->sanitiseName();
	mass_ = new LauParameter(sanitisedName_+"_MASS",mass,0.0,3.0*mass,kTRUE);
	+ mass_->secondStage(kTRUE);
	width_ = new LauParameter(sanitisedName_+"_WIDTH",width,0.0,3.0*width,kTRUE);
	+ width_->secondStage(kTRUE);
	}

	LauResonanceInfo::~LauResonanceInfo()
	{
	delete mass_; mass_ = 0;
	delete width_; width_ = 0;
	}

	/*
	LauResonanceInfo::LauResonanceInfo( const LauResonanceInfo& other ) :
	name_( other.name_ ),
	sanitisedName_( other.sanitisedName_ ),
	mass_( other.mass_->createClone() ),
	width_( other.width_->createClone() ),
	spin_( other.spin_ ),
	charge_( other.charge_ ),
	range_( other.range_ ),
	conjugate_( other.conjugate_ ),
	extraPars_( other.extraPars_ )
	{
	}
	*/

	LauResonanceInfo::LauResonanceInfo( const LauResonanceInfo& other, const TString& newName, const Int_t newCharge ) :
	name_(newName),
	sanitisedName_(""),
	mass_(0),
	width_(0),
	spin_(other.spin_),
	charge_(newCharge),
	range_(other.range_),
	conjugate_(0),
	extraPars_()
	{
	this->sanitiseName();
	mass_ = other.mass_->createClone( sanitisedName_+"_MASS" );
	width_ = other.mass_->createClone( sanitisedName_+"_WIDTH" );
	for ( std::set<LauParameter*>::iterator iter = other.extraPars_.begin(); iter != other.extraPars_.end(); ++iter ) {
	TString parName = (*iter)->name();
	parName.Remove(0, parName.Last('_'));
	parName.Prepend( sanitisedName_ );
	LauParameter* par = (*iter)->createClone( parName );
	extraPars_.insert( par );
	}
	}

	LauResonanceInfo* LauResonanceInfo::createChargeConjugate()
	{
	Int_t newCharge = -charge_;

	TString newName( name_ );
	Ssiz_t index = newName.Index("+");
	if ( index != -1 ) {
	newName.Replace( index, 1, "-" );
	} else {
	index = newName.Index("-");
	if ( index != -1 ) {
	newName.Replace( index, 1, "+" );
	}
	}

	LauResonanceInfo* conjugate = new LauResonanceInfo( *this, newName, newCharge );
	conjugate->conjugate_ = this;
	this->conjugate_ = conjugate;

	return conjugate;
	}

	LauParameter* LauResonanceInfo::getExtraParameter( const TString& parName )
	{
	LauParameter* par(0);
	for ( std::set<LauParameter*>::iterator iter = extraPars_.begin(); iter != extraPars_.end(); ++iter ) {
	if ( (*iter)->name() == parName ) {
	par = (*iter);
	}
	}
	return par;
	}

	void LauResonanceInfo::addExtraParameter( LauParameter* param )
	{
	bool ok = extraPars_.insert( param ).second;
	if ( !ok ) {
	std::cerr << "WARNING in LauResonanceInfo::addExtraParameter : parameter already present, not adding again" << std::endl;
	return;
	}

	if ( conjugate_ != 0 ) {
	conjugate_->addCloneOfExtraParameter(param);
	}
	}

	void LauResonanceInfo::addCloneOfExtraParameter( LauParameter* param )
	{
	TString parName = param->name();
	parName.Remove(0, parName.Last('_'));
	parName.Prepend( sanitisedName_ );

	LauParameter* cloneParam = param->createClone( parName );
	extraPars_.insert( cloneParam );
	}

	ostream& operator<<( ostream& stream, const LauResonanceInfo& infoRecord )
	{
	stream << infoRecord.getName() << ": ";
	stream << "mass = " << infoRecord.getMass()->value() << ", ";
	stream << "width = " << infoRecord.getWidth()->value() << ", ";
	stream << "spin = " << infoRecord.getSpin() << ", ";
	Int_t charge = infoRecord.getCharge();
	if ( charge < 0 ) {
	stream << "charge = " << infoRecord.getCharge() << ", ";
	} else {
	stream << "charge = " << infoRecord.getCharge() << ", ";
	}
	stream << "range = " << infoRecord.getRange();

	return stream;
	}

	void LauResonanceInfo::sanitiseName()
	{
	sanitisedName_ = name_;
	sanitisedName_ = sanitisedName_.ReplaceAll("+","p");
	sanitisedName_ = sanitisedName_.ReplaceAll("-","m");
	sanitisedName_ = sanitisedName_.ReplaceAll("*","st");
	sanitisedName_ = sanitisedName_.ReplaceAll("(","_");
	sanitisedName_ = sanitisedName_.ReplaceAll(")","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("[","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("]","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("<","_");
	sanitisedName_ = sanitisedName_.ReplaceAll(">","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("{","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("}","_");
	sanitisedName_ = sanitisedName_.ReplaceAll(" ","_");
	sanitisedName_ = sanitisedName_.ReplaceAll("$","");
	sanitisedName_ = sanitisedName_.ReplaceAll("%","");
	sanitisedName_ = sanitisedName_.ReplaceAll("&","");
	sanitisedName_ = sanitisedName_.ReplaceAll("/","");
	sanitisedName_ = sanitisedName_.ReplaceAll(":","");
	sanitisedName_ = sanitisedName_.ReplaceAll(";","");
	sanitisedName_ = sanitisedName_.ReplaceAll("=","");
	sanitisedName_ = sanitisedName_.ReplaceAll("\\","");
	sanitisedName_ = sanitisedName_.ReplaceAll("^","");
	sanitisedName_ = sanitisedName_.ReplaceAll("\|","");
	sanitisedName_ = sanitisedName_.ReplaceAll(",","");
	sanitisedName_ = sanitisedName_.ReplaceAll(".","");
	sanitisedName_.Remove(TString::kBoth,'_');
	}
	diff --git a/src/LauSigmaRes.cc b/src/LauSigmaRes.cc
	index 2485a8a..2ce88b5 100644
	--- a/src/LauSigmaRes.cc
	+++ b/src/LauSigmaRes.cc
	@@ -1,288 +1,292 @@

	// Copyright University of Warwick 2004 - 2014.
	// Distributed under the Boost Software License, Version 1.0.
	// (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

	// Authors:
	// Thomas Latham
	// John Back
	// Paul Harrison

	/*! \file LauSigmaRes.cc
	\brief File containing implementation of LauSigmaRes class.
	*/

	#include <iostream>

	#include "LauConstants.hh"
	#include "LauSigmaRes.hh"
	#include "LauResonanceInfo.hh"

	ClassImp(LauSigmaRes)


	LauSigmaRes::LauSigmaRes(LauResonanceInfo* resInfo, const Int_t resPairAmpInt, const LauDaughters* daughters) :
	LauAbsResonance(resInfo, resPairAmpInt, daughters),
	mPiSq4_(4.0*LauConstants::mPiSq),
	sAdler_(LauConstants::mPiSq*0.5),
	b1_(0),
	b2_(0),
	a_(0),
	m0_(0)
	{
	// Initialise various constants
	mPiSq4_ = 4.0*LauConstants::mPiSq;
	sAdler_ = LauConstants::mPiSq0.5; // Adler zero at 0.5(mpi)^2

	// constant factors from BES data
	const Double_t b1Val = 0.5843;
	const Double_t b2Val = 1.6663;
	const Double_t aVal = 1.082;
	const Double_t m0Val = 0.9264;

	const TString& parNameBase = this->getSanitisedName();

	TString b1Name(parNameBase);
	b1Name += "_b1";
	b1_ = resInfo->getExtraParameter( b1Name );
	if ( b1_ == 0 ) {
	b1_ = new LauParameter( b1Name, b1Val, 0.0, 100.0, kTRUE );
	+ b1_->secondStage(kTRUE);
	resInfo->addExtraParameter( b1_ );
	}

	TString b2Name(parNameBase);
	b2Name += "_b2";
	b2_ = resInfo->getExtraParameter( b2Name );
	if ( b2_ == 0 ) {
	b2_ = new LauParameter( b2Name, b2Val, 0.0, 100.0, kTRUE );
	+ b2_->secondStage(kTRUE);
	resInfo->addExtraParameter( b2_ );
	}

	TString aName(parNameBase);
	aName += "_A";
	a_ = resInfo->getExtraParameter( aName );
	if ( a_ == 0 ) {
	a_ = new LauParameter( aName, aVal, 0.0, 10.0, kTRUE );
	+ a_->secondStage(kTRUE);
	resInfo->addExtraParameter( a_ );
	}

	TString m0Name(parNameBase);
	m0Name += "_m0";
	m0_ = resInfo->getExtraParameter( m0Name );
	if ( m0_ == 0 ) {
	m0_ = new LauParameter( m0Name, m0Val, 0.0, 10.0, kTRUE );
	+ m0_->secondStage(kTRUE);
	resInfo->addExtraParameter( m0_ );
	}
	}

	LauSigmaRes::~LauSigmaRes()
	{
	}

	void LauSigmaRes::initialise()
	{
	this->checkDaughterTypes();

	Double_t resSpin = this->getSpin();
	if (resSpin != 0) {
	std::cerr << "WARNING in LauSigmaRes::initialise : Resonance spin is " << resSpin << "." << std::endl;
	std::cerr << " : Sigma amplitude is only for scalers, resetting spin to 0." << std::endl;
	this->changeResonance( -1.0, -1.0, 0 );
	}
	}

	void LauSigmaRes::checkDaughterTypes() const
	{
	// Check that the daughter tracks are the same type. Otherwise issue a warning
	// and set the type to be pion for the Sigma distribution.
	Int_t resPairAmpInt = this->getPairInt();
	if (resPairAmpInt < 1 \|\| resPairAmpInt > 3) {
	std::cerr << "WARNING in LauSigmaRes::checkDaughterTypes : resPairAmpInt = " << resPairAmpInt << " is out of the range [1,2,3]." << std::endl;
	return;
	}

	const TString& nameDaug1 = this->getNameDaug1();
	const TString& nameDaug2 = this->getNameDaug2();
	if (!nameDaug1.CompareTo(nameDaug2, TString::kExact)) {
	// Daughter types agree. Find out if we have pion or kaon system
	if (!nameDaug1.Contains("pi")) {
	std::cerr << "ERROR in LauSigmaRes::checkDaughterTypes : Sigma model is using daughters \""<<nameDaug1<<"\" and \""<<nameDaug2<<"\", which are not pions." << std::endl;
	}
	}
	}

	LauComplex LauSigmaRes::resAmp(Double_t mass, Double_t spinTerm)
	{
	// This function returns the complex dynamical amplitude for a Sigma distribution
	// given the invariant mass and cos(helicity) values.

	// First check that the appropriate daughters are either pi+pi- or K+K-
	// Check that the daughter tracks are the same type. Otherwise issue a warning
	// and set the type to be pion for the Sigma distribution. Returns the
	// integer defined by the enum LauSigmaRes::SigmaPartType.

	Double_t s = mass*mass; // Invariant mass squared combination for the system
	Double_t rho(0.0); // Phase-space factor
	if (s > mPiSq4_) {rho = TMath::Sqrt(1.0 - mPiSq4_/s);}

	const Double_t m0Val = this->getM0Value();
	const Double_t m0Sq = m0Val * m0Val;

	const Double_t aVal = this->getAValue();
	const Double_t b1Val = this->getB1Value();
	const Double_t b2Val = this->getB2Value();

	Double_t f = b2Val*s + b1Val; // f(s) function
	Double_t numerator = s - sAdler_;
	Double_t denom = m0Sq - sAdler_;
	Double_t gamma(0.0);
	if (TMath::Abs(denom) > 1e-10 && TMath::Abs(aVal) > 1e-10) {
	// Decay width of the system
	gamma = rho(numerator/denom)f*TMath::Exp(-(s - m0Sq)/aVal);
	}

	// Now form the complex amplitude - use relativistic BW form (without barrier factors)
	// Note that the M factor in the denominator is not the "pole" at ~500 MeV, but is
	// m0_ = 0.9264, the mass when the phase shift goes through 90 degrees.

	Double_t dMSq = m0Sq - s;
	Double_t widthTerm = gamma*m0Val;
	LauComplex resAmplitude(dMSq, widthTerm);

	Double_t denomFactor = dMSqdMSq + widthTermwidthTerm;

	Double_t invDenomFactor = 0.0;
	if (denomFactor > 1e-10) {invDenomFactor = 1.0/denomFactor;}

	resAmplitude.rescale(spinTerm*invDenomFactor);

	return resAmplitude;

	}

	const std::vector<LauParameter*>& LauSigmaRes::getFloatingParameters()
	{
	this->clearFloatingParameters();

	if ( ! this->fixB1Value() ) {
	this->addFloatingParameter( b1_ );
	}

	if ( ! this->fixB2Value() ) {
	this->addFloatingParameter( b2_ );
	}

	if ( ! this->fixAValue() ) {
	this->addFloatingParameter( a_ );
	}

	if ( ! this->fixM0Value() ) {
	this->addFloatingParameter( m0_ );
	}

	return this->getParameters();
	}

	void LauSigmaRes::setResonanceParameter(const TString& name, const Double_t value)
	{
	// Set various parameters for the lineshape
	if (name == "b1") {
	this->setB1Value(value);
	std::cout << "INFO in LauSigmaRes::setResonanceParameter : Setting parameter b1 = " << this->getB1Value() << std::endl;
	}
	else if (name == "b2") {
	this->setB2Value(value);
	std::cout << "INFO in LauSigmaRes::setResonanceParameter : Setting parameter b2 = " << this->getB2Value() << std::endl;
	}
	else if (name == "A") {
	this->setAValue(value);
	std::cout << "INFO in LauSigmaRes::setResonanceParameter : Setting parameter A = " << this->getAValue() << std::endl;
	}
	else if (name == "m0") {
	this->setM0Value(value);
	std::cout << "INFO in LauSigmaRes::setResonanceParameter : Setting parameter m0 = " << this->getM0Value() << std::endl;
	}
	else {
	std::cerr << "WARNING in LauSigmaRes::setResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	void LauSigmaRes::floatResonanceParameter(const TString& name)
	{
	if (name == "b1") {
	if ( b1_->fixed() ) {
	b1_->fixed( kFALSE );
	this->addFloatingParameter( b1_ );
	} else {
	std::cerr << "WARNING in LauSigmaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "b2") {
	if ( b2_->fixed() ) {
	b2_->fixed( kFALSE );
	this->addFloatingParameter( b2_ );
	} else {
	std::cerr << "WARNING in LauSigmaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "A") {
	if ( a_->fixed() ) {
	a_->fixed( kFALSE );
	this->addFloatingParameter( a_ );
	} else {
	std::cerr << "WARNING in LauSigmaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else if (name == "m0") {
	if ( m0_->fixed() ) {
	m0_->fixed( kFALSE );
	this->addFloatingParameter( m0_ );
	} else {
	std::cerr << "WARNING in LauSigmaRes::floatResonanceParameter: Parameter already floating. No parameter changes made." << std::endl;
	}
	} else {
	std::cerr << "WARNING in LauSigmaRes::fixResonanceParameter: Parameter name not reconised. No parameter changes made." << std::endl;
	}
	}

	LauParameter* LauSigmaRes::getResonanceParameter(const TString& name)
	{
	if (name == "b1") {
	return b1_;
	} else if (name == "b2") {
	return b2_;
	} else if (name == "A") {
	return a_;
	} else if (name == "m0") {
	return m0_;
	} else {
	std::cerr << "WARNING in LauSigmaRes::getResonanceParameter: Parameter name not reconised." << std::endl;
	return 0;
	}
	}

	void LauSigmaRes::setB1Value(const Double_t b1)
	{
	b1_->value( b1 );
	b1_->genValue( b1 );
	b1_->initValue( b1 );
	}

	void LauSigmaRes::setB2Value(const Double_t b2)
	{
	b2_->value( b2 );
	b2_->genValue( b2 );
	b2_->initValue( b2 );
	}

	void LauSigmaRes::setAValue(const Double_t A)
	{
	a_->value( A );
	a_->genValue( A );
	a_->initValue( A );
	}

	void LauSigmaRes::setM0Value(const Double_t m0)
	{
	m0_->value( m0 );
	m0_->genValue( m0 );
	m0_->initValue( m0 );
	}

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 3:41 PM (1 d, 21 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805000
Default Alt Text: (83 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions