Page MenuHomeHEPForge
Files F7877532

No OneTemporary
Actions

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.0*nSigEvents, 2.0*nSigEvents, 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(-alpha*mass*mass);
} 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(-alpha*s) + TMath::Exp(-alpha*t);
} 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.5*TMath::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 = realPart*realPart + imagPart*imagPart;
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 != "K*0_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 += g2Val*resMass*0.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 += g2Val*resMass*(0.5*TMath::Sqrt(mSumSq2_/s - 1.0) + 0.5*TMath::Sqrt(mSumSq3_/s - 1.0));
}
} else {
// Continue analytically below higher channel thresholds
// This contributes to the real part of the amplitude denominator
dMSq += g1Val*resMass*2.0*TMath::Sqrt(mSumSq1_/s - 1.0)/3.0;
}
}
Double_t massFactor = resMass;
if(useAdlerTerm_) massFactor = ( s - sA_ ) / ( resMassSq - sA_ );
const Double_t width1 = g1Val*rho1*massFactor;
const Double_t width2 = g2Val*rho2*massFactor;
const Double_t widthTerm = width1 + width2;
LauComplex resAmplitude(dMSq, widthTerm);
const Double_t denomFactor = dMSq*dMSq + widthTerm*widthTerm;
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 = dMSq*dMSq + widthTerm*widthTerm;
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.0*resMass_);
}
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 = resWidth*qRatio*(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_)/(massSqTerm*massSqTerm + resMassSq_*totWidth*totWidth));
// Calculate the phase shift term
const Double_t rVal = this->getEffectiveRange();
const Double_t aVal = this->getScatteringLength();
const Double_t tandeltaB = (2.0*aVal*q)/(2.0 + aVal*rVal*q*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 + (rVal*q*q)/2.0;
// Compute the complex amplitude
bkgAmplitude = LauComplex(qcotdeltaB, q);
// Scale by the numerator and denominator factors
bkgAmplitude.rescale(mass/(qcotdeltaB*qcotdeltaB + q*q));
}
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.0*resMass_);
}
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 = resWidth*qRatio*(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_)/(massSqTerm*massSqTerm + resMassSq_*totWidth*totWidth));
// Calculate the phase shift term
const Double_t rVal = this->getEffectiveRange();
const Double_t aVal = this->getScatteringLength();
const Double_t tandeltaB = (2.0*aVal*q)/(2.0 + aVal*rVal*q*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 + (rVal*q*q)/2.0;
// Compute the complex amplitude
bkgAmplitude = LauComplex(qcotdeltaB, q);
// Scale by the numerator and denominator factors
bkgAmplitude.rescale(mass/(qcotdeltaB*qcotdeltaB + q*q));
//------------------
// 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() * m13Sq*m13Sq * 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::mPiSq*0.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 = dMSq*dMSq + widthTerm*widthTerm;
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, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805000
Default Alt Text
(83 KB)

Event Timeline