LauFlatteRes.cc

Go to the documentation of this file.

// 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

#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 - Phys. Lett. B 607, 243 (2005)
                // resMass should be 0.965 +/- 0.008 +/- 0.006 GeV/c^2
                g1Val = 0.165;       // +/- 0.010 +/- 0.015 GeV^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()
{
}

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 );
}