LauAbsModIndPartWave.cc

Go to the documentation of this file.

// Copyright University of Warwick 2004 - 2013.
// 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 <cstdlib>
#include <iostream>

#include "LauConstants.hh"
#include "LauKinematics.hh"
#include "LauAbsModIndPartWave.hh"
#include "LauResonanceInfo.hh"

ClassImp(LauAbsModIndPartWave)


LauAbsModIndPartWave::LauAbsModIndPartWave(LauResonanceInfo* resInfo, Int_t resPairAmpInt, const LauDaughters* daughters) :
        LauAbsResonance(resInfo, resPairAmpInt, daughters),
        nKnots_(0),
        spline1_(0),
        spline2_(0),
        type1_(Lau1DCubicSpline::StandardSpline),
        type2_(Lau1DCubicSpline::StandardSpline),
        leftBound1_(Lau1DCubicSpline::NotAKnot),
        rightBound1_(Lau1DCubicSpline::NotAKnot),
        leftBound2_(Lau1DCubicSpline::NotAKnot),
        rightBound2_(Lau1DCubicSpline::NotAKnot),
        leftGrad1_(0.),
        rightGrad1_(0.),
        leftGrad2_(0.),
        rightGrad2_(0.),
        secondStage_(kFALSE)
{
}

LauAbsModIndPartWave::~LauAbsModIndPartWave()
{
        delete spline1_;
        delete spline2_;
}

void LauAbsModIndPartWave::floatKnotsSecondStage(const Bool_t secondStage)
{
        secondStage_ = secondStage;

        // if the parameters have not yet been created we can now just return
        if ( amp1Pars_.size() != nKnots_ ) {
                return;
        }

        // otherwise we need to toggle their second-stage parameter
        for ( UInt_t i(0); i < nKnots_; ++i ) {
                amp1Pars_[i]->secondStage(secondStage_);
                amp2Pars_[i]->secondStage(secondStage_);
        }
}

std::set<Double_t> LauAbsModIndPartWave::checkKnots(const std::set<Double_t>& masses)
{
        std::set<Double_t> knots = masses;

        const std::set<Double_t>::const_iterator first = knots.begin();
        const std::set<Double_t>::const_reverse_iterator last = knots.rbegin();

        const Double_t lower_limit = this->getMassDaug1() + this->getMassDaug2();
        const Double_t upper_limit = this->getMassParent() - this->getMassBachelor();

        // check whether we have been given knots at unphysical masses
        if ( *first < lower_limit ) {
                std::cerr << "WARNING in LauAbsModIndPartWave::checkKnots : Knot found at mass " << *first << " is below the lower kinematic limit." << std::endl;
                std::cerr << "                                            : Lower kinematic limit is at mass " << lower_limit << std::endl;
                std::cerr << "                                            : Aborting definition of knot positions." << std::endl;
                knots.clear();
                return knots;
        }
        if ( *last > upper_limit ) {
                std::cerr << "WARNING in LauAbsModIndPartWave::checkKnots : Knot found at mass " << *last << " is above the upper kinematic limit." << std::endl;
                std::cerr << "                                            : Upper kinematic limit is at mass " << upper_limit << std::endl;
                std::cerr << "                                            : Aborting definition of knot positions." << std::endl;
                knots.clear();
                return knots;
        }

        // check if we have knots at each extreme - if not, add them in
        if ( (*first) != lower_limit ) {
                knots.insert( lower_limit );
        }
        if ( (*last)  != upper_limit ) {
                knots.insert( upper_limit );
        }

        return knots;
}

void LauAbsModIndPartWave::defineKnots(const std::set<Double_t>& masses)
{
        if ( ! masses_.empty() ) {
                std::cerr << "WARNING in LauAbsModIndPartWave::defineKnots : Knot positions have already been defined, not making any changes." << std::endl;
                return;
        }

        const std::set<Double_t> knots = this->checkKnots( masses );

        nKnots_ = knots.size();
        if ( nKnots_ == 0 ) {
                return;
        }

        masses_.reserve(nKnots_);
        amp1Vals_.reserve(nKnots_);
        amp2Vals_.reserve(nKnots_);
        amp1Pars_.reserve(nKnots_);
        amp2Pars_.reserve(nKnots_);

        UInt_t counter(0);
        for ( std::set<Double_t>::const_iterator iter = knots.begin(); iter != knots.end(); ++iter ) {
                masses_.push_back( *iter );
                amp1Vals_.push_back(1.0);
                amp2Vals_.push_back(1.0);
                this->createAmpParameters(counter);
                ++counter;
        }

        for ( std::vector<Double_t>::const_iterator iter = masses_.begin(); iter != masses_.end(); ++iter ) {
                std::cout << "INFO in LauAbsModIndPartWave::defineKnots : Knot added to resonance " << this->getResonanceName() << " at mass " << *iter << std::endl;
        }
}

void LauAbsModIndPartWave::initialise()
{
        if ( spline1_ != 0 ) {
                delete spline1_;
                spline1_ = 0;
        }
        if ( spline2_ != 0 ) {
                delete spline2_;
                spline2_ = 0;
        }

        for ( UInt_t i(0); i < nKnots_; ++i ) {
                amp1Vals_[i] = amp1Pars_[i]->unblindValue();
                amp2Vals_[i] = amp2Pars_[i]->unblindValue();
        }

        spline1_ = new Lau1DCubicSpline(masses_, amp1Vals_, type1_, leftBound1_, rightBound1_, leftGrad1_, rightGrad1_);
        spline2_ = new Lau1DCubicSpline(masses_, amp2Vals_, type2_, leftBound2_, rightBound2_, leftGrad2_, rightGrad2_);
}

LauComplex LauAbsModIndPartWave::resAmp(Double_t mass, Double_t spinTerm)
{
        amp_.zero();

        Bool_t paramChanged1(kFALSE), paramChanged2(kFALSE);

        for ( UInt_t i(0); i < nKnots_; ++i ) {
                if ( !amp1Pars_[i]->fixed() && amp1Pars_[i]->unblindValue() != amp1Vals_[i] ) {
                        paramChanged1 = kTRUE;
                        amp1Vals_[i] = amp1Pars_[i]->unblindValue();
                }
                if ( !amp2Pars_[i]->fixed() && amp2Pars_[i]->unblindValue() != amp2Vals_[i] ) {
                        paramChanged2 = kTRUE;
                        amp2Vals_[i] = amp2Pars_[i]->unblindValue();
                }
        }

        if ( spline1_ == 0 ||  spline2_ == 0) {
                std::cerr << "ERROR in LauAbsModIndPartWave::resAmp : One of the splines is null" << std::endl;
                return amp_;
        }

        if ( paramChanged1 ) {
                spline1_->updateYValues(amp1Vals_);
        }

        if ( paramChanged2 ) {
                spline2_->updateYValues(amp2Vals_);
        }

        this->evaluateAmplitude( mass );

        amp_.rescale( spinTerm );

        return amp_;
}

void LauAbsModIndPartWave::setSplineType(Lau1DCubicSpline::LauSplineType type1, Lau1DCubicSpline::LauSplineType type2)
{
        type1_ = type1;
        type2_ = type2;
}

void LauAbsModIndPartWave::setSplineBoundaryConditions(Lau1DCubicSpline::LauSplineBoundaryType leftBound1, 
                                                       Lau1DCubicSpline::LauSplineBoundaryType rightBound1,
                                                       Lau1DCubicSpline::LauSplineBoundaryType leftBound2,
                                                       Lau1DCubicSpline::LauSplineBoundaryType rightBound2,
                                                       Double_t leftGrad1,
                                                       Double_t rightGrad1,
                                                       Double_t leftGrad2,
                                                       Double_t rightGrad2)
{
        leftBound1_  = leftBound1;
        rightBound1_ = rightBound1;
        leftBound2_  = leftBound2;
        rightBound2_ = rightBound2;
        leftGrad1_   = leftGrad1;
        rightGrad1_  = rightGrad1;
        leftGrad2_   = leftGrad2;
        rightGrad2_  = rightGrad2;
}

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

        for ( UInt_t i(0); i < nKnots_; ++i ) {
                if ( !amp1Pars_[i]->fixed() ) {
                        this->addFloatingParameter( amp1Pars_[i] );
                }
                if ( !amp2Pars_[i]->fixed() ) {
                        this->addFloatingParameter( amp2Pars_[i] );
                }
        }

        return this->getParameters();
}