LauFlatteRes.cc

Go to the documentation of this file.

/*
Copyright 2004 University of Warwick

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
*/

/*
Laura++ package authors:
John Back
Paul Harrison
Thomas Latham
*/

#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),
        absorbM0_(kFALSE)
{
        Double_t resMass(0.0);
        Double_t g1Val(0.0);
        Double_t g2Val(0.0);

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

        if ( resName == "f_0(980)" ) {

                // Set the mass threshold values
                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);

                // Set values of mass and coupling constants from:
                // Phys. Lett. B 607, 243 (2005)
                resMass = 0.965;      // 0.965 +/- 0.008 +/- 0.006 GeV/c^2
                g1Val   = 0.165;      // 0.165 +/- 0.010 +/- 0.015 GeV^2
                g2Val   = 4.21*g1Val; // 4.21  +/- 0.25  +/- 0.21

                useAdlerTerm_ = kFALSE;
                sA_ = 0.0;

                absorbM0_ = kTRUE;

        } else if ( resName == "K*0_0(1430)" ) {

                // Set the mass threshold values
                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);

                // Set values of mass and coupling constants from:
                // Phys. Lett. B 572, 1 (2003)
                resMass = 1.513; // GeV/c^2
                g1Val   = 0.304; // GeV
                g2Val   = 0.380; // GeV

                useAdlerTerm_ = kTRUE;
                sA_ = 0.234;

                absorbM0_ = kFALSE;

        } else if ( resName == "K*+_0(1430)" || resName == "K*-_0(1430)" ) {

                // Set the mass threshold values
                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);

                // Set values of mass and coupling constants from:
                // Phys. Lett. B 572, 1 (2003)
                resMass = 1.513; // GeV/c^2
                g1Val   = 0.304; // GeV
                g2Val   = 0.380; // GeV

                useAdlerTerm_ = kTRUE;
                sA_ = 0.234;

                absorbM0_ = kFALSE;

        } else if ( resName == "a0_0(980)" ) {

                // Set the mass threshold values
                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);

                // Set values of mass and coupling constants from:
                // Phys. Rev. D 57, 3860 (1998)
                resMass = 0.982;       // 0.982 +/- 0.003 GeV/c^2
                g1Val   = 0.324*0.324; // 0.324 +/- 0.015 GeV     (NB this value needs to be squared since the paper uses g_1^2 and g_2^2)
                g2Val   = 1.03*g1Val;  // 1.03  +/- 0.14          (NB this is indeed the ratio of what the paper refers to as g_2^2 and g_1^2, so can be used unchanged here)

                useAdlerTerm_ = kFALSE;
                sA_ = 0.0;

                absorbM0_ = kTRUE;

        } else if ( resName == "a+_0(980)" || resName == "a-_0(980)" ) {

                // Set the mass threshold values
                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);

                // Set values of mass and coupling constants from:
                // Phys. Rev. D 57, 3860 (1998)
                resMass = 0.982;       // 0.982 +/- 0.003 GeV/c^2
                g1Val   = 0.324*0.324; // 0.324 +/- 0.015 GeV     (NB this value needs to be squared since the paper uses g_1^2 and g_2^2)
                g2Val   = 1.03*g1Val;  // 1.03  +/- 0.14          (NB this is indeed the ratio of what the paper refers to as g_2^2 and g_1^2, so can be used unchanged here)

                useAdlerTerm_ = kFALSE;
                sA_ = 0.0;

                absorbM0_ = kTRUE;

        }

        const TString couplingUnits = (absorbM0_) ? "GeV^2" : "GeV";
        std::cout << "INFO in LauFlatteRes::LauFlatteRes : Setting default parameters for " << resName << ":\n";
        std::cout << "                                   : mass = " << resMass << " GeV/c^2\n";
        std::cout << "                                   : g1   = " << g1Val << " " << couplingUnits << "\n";
        std::cout << "                                   : g2   = " << g2Val << " " << couplingUnits << "\n";
        if ( absorbM0_ ) {
                std::cout << "                                   : Will absorb m0 into couplings\n";
        } else {
                std::cout << "                                   : Will not absorb m0 into couplings\n";
        }
        if ( useAdlerTerm_ ) {
                std::cout << "                                   : Will use Adler zero term\n";
                std::cout << "                                   : sA   = " << sA_ << " GeV^2/c^4\n";
        } else {
                std::cout << "                                   : Will not use Adler zero term\n";
        }

        // Set the mass value
        LauParameter* massPar = this->getMassPar();
        if ( massPar ) {
                massPar->valueAndRange(resMass,0.0,3.0*resMass);
                massPar->initValue(resMass);
                massPar->genValue(resMass);
        } else {
                std::cerr << "ERROR in LauFlatteRes::LauFlatteRes : Unable to retrieve mass parameter" << std::endl;
        }

        // Create the parameters for the coupling constants
        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()
{
}

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 = 1.0;
        if ( ! absorbM0_ ) {
                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 );
}