LauParameter.cc

Go to the documentation of this file.

// Copyright University of Warwick 2006 - 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 <iostream>
#include <map>
using std::cout;
using std::cerr;
using std::endl;
using std::map;

#include "TRandom.h"

#include "LauParameter.hh"
#include "LauRandom.hh"

ClassImp(LauParameter)


LauParameter::LauParameter() :
        name_(""),
        value_(0.0),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(0.0),
        initValue_(0.0),
        minValue_(0.0),
        maxValue_(0.0),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
}

LauParameter::LauParameter(const TString& parName) :
        name_(parName),
        value_(0.0),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(0.0),
        initValue_(0.0),
        minValue_(0.0),
        maxValue_(0.0),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
}

LauParameter::LauParameter(Double_t parValue) :
        name_(""),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(parValue-1e-6),
        maxValue_(parValue+1e-6),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
}

LauParameter::LauParameter(const TString& parName, Double_t parValue) :
        name_(parName),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(parValue-1e-6),
        maxValue_(parValue+1e-6),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
}

LauParameter::LauParameter(Double_t parValue, Double_t min, Double_t max) :
        name_(""),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(Double_t parValue, Double_t parError, Double_t min, Double_t max) :
        name_(""),
        value_(parValue),
        error_(parError),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(Double_t parValue, Double_t min, Double_t max, Bool_t parFixed) :
        name_(""),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(parFixed),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t min, Double_t max) :
        name_(parName),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t min, Double_t max, Bool_t parFixed) :
        name_(parName),
        value_(parValue),
        error_(0.0),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(parFixed),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(const TString& parName, Double_t parValue, Double_t parError, Double_t min, Double_t max) :
        name_(parName),
        value_(parValue),
        error_(parError),
        negError_(0.0),
        posError_(0.0),
        genValue_(parValue),
        initValue_(parValue),
        minValue_(min),
        maxValue_(max),
        fixed_(kTRUE),
        firstStage_(kFALSE),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0)
{
        this->checkRange();
}

LauParameter::LauParameter(const LauParameter& rhs) : TObject(rhs), LauAbsRValue(rhs),
        name_(rhs.name_),
        value_(rhs.value_),
        error_(rhs.error_),
        negError_(rhs.negError_),
        posError_(rhs.posError_),
        genValue_(rhs.genValue_),
        initValue_(rhs.initValue_),
        minValue_(rhs.minValue_),
        maxValue_(rhs.maxValue_),
        fixed_(rhs.fixed_),
        firstStage_(rhs.firstStage_),
        secondStage_(rhs.secondStage_),
        gaussConstraint_(rhs.gaussConstraint_),
        constraintMean_(rhs.constraintMean_),
        constraintWidth_(rhs.constraintWidth_),
        gcc_(rhs.gcc_),
        bias_(rhs.bias_),
        pull_(rhs.pull_),
        clone_(rhs.clone_),
        parent_(rhs.parent_),
        clones_(rhs.clones_)
{
}

LauParameter& LauParameter::operator=(const LauParameter& rhs)
{
        if (&rhs != this) {
                TObject::operator=(rhs);
                LauAbsRValue::operator=(rhs);
                name_ = rhs.name_;
                value_ = rhs.value_;
                error_ = rhs.error_;
                negError_ = rhs.negError_;
                posError_ = rhs.posError_;
                genValue_ = rhs.genValue_;
                initValue_ = rhs.initValue_;
                minValue_ = rhs.minValue_;
                maxValue_ = rhs.maxValue_;
                fixed_ = rhs.fixed_;
                firstStage_ = rhs.firstStage_;
                secondStage_ = rhs.secondStage_;
                gaussConstraint_ = rhs.gaussConstraint_;
                constraintMean_ = rhs.constraintMean_;
                constraintWidth_ = rhs.constraintWidth_;
                gcc_ = rhs.gcc_;
                bias_ = rhs.bias_;
                pull_ = rhs.pull_;
                clone_ = rhs.clone_;
                parent_ = rhs.parent_;
                clones_ = rhs.clones_;
        }
        return *this;
}

std::vector<LauParameter*> LauParameter::getPars() 
{
        std::vector<LauParameter*> list;
        list.push_back(this);
        return list;
}

void LauParameter::value(Double_t newValue)
{
        if (this->clone()) {
                parent_->value(newValue);
        } else {
                this->checkRange(newValue,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
}

void LauParameter::error(Double_t newError)
{
        if (this->clone()) {
                parent_->error(newError);
        } else {
                error_ = TMath::Abs(newError);
                this->updateClones(kFALSE);
        }
}

void LauParameter::negError(Double_t newNegError)
{
        if (this->clone()) {
                parent_->negError(newNegError);
        } else {
                negError_ = TMath::Abs(newNegError);
                this->updateClones(kFALSE);
        }
}

void LauParameter::posError(Double_t newPosError)
{
        if (this->clone()) {
                parent_->posError(newPosError);
        } else {
                posError_ = TMath::Abs(newPosError);
                this->updateClones(kFALSE);
        }
}

void LauParameter::errors(Double_t newError, Double_t newNegError, Double_t newPosError)
{
        if (this->clone()) {
                parent_->errors(newError,newNegError,newPosError);
        } else {
                error_ = TMath::Abs(newError);
                negError_ = TMath::Abs(newNegError);
                posError_ = TMath::Abs(newPosError);
                this->updateClones(kFALSE);
        }
}

void LauParameter::valueAndErrors(Double_t newValue, Double_t newError, Double_t newNegError, Double_t newPosError)
{
        if (this->clone()) {
                parent_->valueAndErrors(newValue,newError,newNegError,newPosError);
        } else {
                this->checkRange(newValue,this->minValue(),this->maxValue());
                error_ = TMath::Abs(newError);
                negError_ = TMath::Abs(newNegError);
                posError_ = TMath::Abs(newPosError);
                this->updateClones(kFALSE);
        }
}

void LauParameter::globalCorrelationCoeff(Double_t newGCCValue)
{
        if (this->clone()) {
                parent_->globalCorrelationCoeff(newGCCValue);
        } else {
                gcc_ = newGCCValue;
                this->updateClones(kFALSE);
        }
}

void LauParameter::genValue(Double_t newGenValue)
{
        if (this->clone()) {
                parent_->genValue(newGenValue);
        } else {
                genValue_ = newGenValue;
                this->updateClones(kFALSE);
        }
}

void LauParameter::initValue(Double_t newInitValue)
{
        if (this->clone()) {
                parent_->initValue(newInitValue);
        } else {
                initValue_ = newInitValue;
                this->updateClones(kFALSE);
        }
}

void LauParameter::minValue(Double_t newMinValue)
{
        if (this->clone()) {
                parent_->minValue(newMinValue);
        } else {
                this->checkRange(this->value(),newMinValue,this->maxValue());
                this->updateClones(kFALSE);
        }
}

void LauParameter::maxValue(Double_t newMaxValue)
{
        if (this->clone()) {
                parent_->maxValue(newMaxValue);
        } else {
                this->checkRange(this->value(),this->minValue(),newMaxValue);
                this->updateClones(kFALSE);
        }
}

void LauParameter::range(Double_t newMinValue, Double_t newMaxValue)
{
        if (this->clone()) {
                parent_->range(newMinValue,newMaxValue);
        } else {
                this->checkRange(this->value(),newMinValue,newMaxValue);
                this->updateClones(kFALSE);
        }
}

void LauParameter::valueAndRange(Double_t newValue, Double_t newMinValue, Double_t newMaxValue)
{
        if (this->clone()) {
                parent_->valueAndRange(newValue,newMinValue,newMaxValue);
        } else {
                this->checkRange(newValue,newMinValue,newMaxValue);
                this->updateClones(kFALSE);
        }
}

void LauParameter::name(const TString& newName)
{
        // no need to update clones here
        // clones are allowed to have different names
        name_ = newName;
}

void LauParameter::fixed(Bool_t parFixed)
{
        if (this->clone()) {
                parent_->fixed(parFixed);
        } else {
                fixed_ = parFixed;
                this->updateClones(kFALSE);
        }
}

void LauParameter::firstStage(Bool_t firstStagePar)
{
        if (this->clone()) {
                parent_->firstStage(firstStagePar);
        } else {
                firstStage_ = firstStagePar;
                this->updateClones(kFALSE);
        }
}

void LauParameter::secondStage(Bool_t secondStagePar)
{
        if (this->clone()) {
                parent_->secondStage(secondStagePar);
        } else {
                secondStage_ = secondStagePar;
                this->updateClones(kFALSE);
        }
}

void LauParameter::addGaussianConstraint(Double_t newGaussMean, Double_t newGaussWidth)
{
        if (this->clone()) {
                parent_->addGaussianConstraint(newGaussMean,newGaussWidth);
        } else {
                gaussConstraint_ = kTRUE;
                constraintMean_ = newGaussMean;
                constraintWidth_ = newGaussWidth;
                this->updateClones(kFALSE);
        }
}

void LauParameter::removeGaussianConstraint()
{
        if (this->clone()) {
                parent_->removeGaussianConstraint();
        } else {
                gaussConstraint_ = kFALSE;
                this->updateClones(kFALSE);
        }
}

LauParameter& LauParameter::operator = (Double_t val)
{
        if (this->clone()) {
                (*parent_) = val;
        } else {
                this->checkRange(val,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
        return *this;
}

LauParameter& LauParameter::operator += (Double_t val)
{
        if (this->clone()) {
                (*parent_) += val;
        } else {
                this->checkRange(this->value()+val,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
        return *this;
}

LauParameter& LauParameter::operator -= (Double_t val)
{
        if (this->clone()) {
                (*parent_) -= val;
        } else {
                this->checkRange(this->value()-val,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
        return *this;
}

LauParameter& LauParameter::operator *= (Double_t val)
{
        if (this->clone()) {
                (*parent_) *= val;
        } else {
                this->checkRange(this->value()*val,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
        return *this;
}

LauParameter& LauParameter::operator /= (Double_t val)
{
        if (this->clone()) {
                (*parent_) /= val;
        } else {
                this->checkRange(this->value()/val,this->minValue(),this->maxValue());
                this->updateClones(kTRUE);
        }
        return *this;
}

void LauParameter::updatePull()
{
        if (this->clone()) {
                parent_->updatePull();
                return;
        }

        // calculate the bias
        bias_ = value_ - genValue_;

        // if we have errors calculated then calculate
        // the pull using the best error available
        if ((bias_ > 0.0) && (negError_ > 1e-10)) {
                pull_ = bias_ / negError_;
        } else if ((bias_ < 0.0) && (posError_ > 1e-10)) {
                pull_ = bias_ / posError_;
        } else if (error_ > 1e-10) {
                pull_ = bias_ / error_;
        } else {
                pull_ = 0.0;
        }

        this->updateClones(kFALSE);
}

void LauParameter::checkRange(Double_t val, Double_t minVal, Double_t maxVal)
{
        // first check that min is less than max (or they are the same - this is allowed)
        if (minVal > maxVal) {
                cerr<<"ERROR in LauParameter::checkRange : minValue: "<<minVal<<" greater than maxValue: "<<maxVal<<endl;
                if (minValue_ > maxValue_) {
                        minValue_ = maxValue_;
                        cerr<<"                                  : Setting both to "<<maxValue_<<"."<<endl;
                } else {
                        cerr<<"                                  : Not changing anything."<<endl;
                }
                return;
        }

        minValue_ = minVal;
        maxValue_ = maxVal;

        // now check that the value is still within the range
        if ((val < minVal) || (val > maxVal)) {
                if (name_ != "") {
                        cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"] for parameter \""<<name_<<"\"."<<endl;
                } else {
                        cerr<<"ERROR in LauParameter::checkRange : value: "<<val<<" not in allowed range ["<<minVal<<" , "<<maxVal<<"]."<<endl;
                }
                if (val < minVal) {
                        cerr<<"                                  : Setting value to minValue: "<<minVal<<endl;
                        val = minVal;
                } else {
                        cerr<<"                                  : Setting value to maxValue: "<<maxVal<<endl;
                        val = maxVal;
                }
        }
        value_ = val;
}

LauParameter* LauParameter::createClone(Double_t constFactor)
{
        // if we're a clone we mustn't be cloned ourselves
        // but instead return another clone of our parent
        // (this is so that the parent knows of all its clones)
        if (this->clone()) {
                LauParameter* clonePar = parent_->createClone(constFactor);
                clonePar->name(this->name());
                return clonePar;
        }

        // clone ourselves using the copy-constructor
        LauParameter* clonePar = new LauParameter(*this);
        (*clonePar) *= constFactor;
        clonePar->wipeClones();
        clonePar->clone(this);
        clones_.insert( std::make_pair( clonePar, constFactor ) );
        return clonePar;
}

LauParameter* LauParameter::createClone(const TString& newName, Double_t constFactor)
{
        // self message to create the clone
        LauParameter* clonePar = this->createClone(constFactor);
        // set the new name
        clonePar->name(newName);
        // and return
        return clonePar;
}

void LauParameter::updateClones(Bool_t justValue)
{
        // if we don't have any clones then there's nothing to do
        if ( clones_.empty() ) {
                return;
        }

        // we have to set the values directly rather than using member functions because otherwise we'd get into an infinite loop
        if (justValue) {
                for (map<LauParameter*,Double_t>::iterator iter = clones_.begin(); iter != clones_.end(); ++iter) {
                        LauParameter* clonePar = iter->first;
                        Double_t constFactor = iter->second;
                        clonePar->value_ = constFactor*this->value();
                }
        } else {
                for (map<LauParameter*,Double_t>::iterator iter = clones_.begin(); iter != clones_.end(); ++iter) {
                        LauParameter* clonePar = iter->first;
                        Double_t constFactor = iter->second;
                        clonePar->value_ = constFactor*this->value();
                        clonePar->error_ = constFactor*this->error();
                        clonePar->negError_ = constFactor*this->negError();
                        clonePar->posError_ = constFactor*this->posError();
                        clonePar->genValue_ = constFactor*this->genValue();
                        clonePar->initValue_ = constFactor*this->initValue();
                        clonePar->minValue_ = constFactor*this->minValue();
                        clonePar->maxValue_ = constFactor*this->maxValue();
                        clonePar->fixed_ = this->fixed();
                        clonePar->firstStage_ = this->firstStage();
                        clonePar->secondStage_ = this->secondStage();
                        clonePar->gaussConstraint_ = this->gaussConstraint();
                        clonePar->constraintMean_ = this->constraintMean();
                        clonePar->constraintWidth_ = this->constraintWidth();
                        clonePar->gcc_ = this->globalCorrelationCoeff();
                        clonePar->bias_ = this->bias();
                        clonePar->pull_ = this->pull();
                }
        }
}

void LauParameter::randomiseValue()
{
        this->randomiseValue(this->minValue(), this->maxValue());
}

void LauParameter::randomiseValue(Double_t minVal, Double_t maxVal)
{
        // if we're fixed then do nothing
        if (this->fixed()) {
                return;
        }

        // check supplied values are sensible
        if (maxVal < minVal) {
                cerr<<"ERROR in LauParameter::randomiseValue : Supplied maximum value smaller than minimum value."<<endl;
                return;
        }
        if (maxVal > this->maxValue()) {
                maxVal = this->maxValue();
        }
        if (minVal < this->minValue()) {
                minVal = this->minValue();
        }

        // use the zero-seed random number generator to get values that are
        // uniformly distributed over the given range
        Double_t randNo = LauRandom::zeroSeedRandom()->Rndm();
        Double_t val = randNo*(maxVal - minVal) + minVal;
        this->initValue(val);
}

// various mathematical operators (non-member functions)
Double_t operator + (const LauParameter& lhs, Double_t rhs)
{
        return (lhs.value() + rhs);
}
Double_t operator + (Double_t lhs, const LauParameter& rhs)
{
        return (lhs + rhs.value());
}
Double_t operator + (const LauParameter& lhs, const LauParameter& rhs)
{
        return (lhs.value() + rhs.value());
}
Double_t operator - (const LauParameter& lhs, Double_t rhs)
{
        return (lhs.value() - rhs);
}
Double_t operator - (Double_t lhs, const LauParameter& rhs)
{
        return (lhs - rhs.value());
}
Double_t operator - (const LauParameter& lhs, const LauParameter& rhs)
{
        return (lhs.value() - rhs.value());
}
Double_t operator * (const LauParameter& lhs, Double_t rhs)
{
        return (lhs.value() * rhs);
}
Double_t operator * (Double_t lhs, const LauParameter& rhs)
{
        return (lhs * rhs.value());
}
Double_t operator * (const LauParameter& lhs, const LauParameter& rhs)
{
        return (lhs.value() * rhs.value());
}
Double_t operator / (const LauParameter& lhs, Double_t rhs)
{
        return (lhs.value() / rhs);
}
Double_t operator / (Double_t lhs, const LauParameter& rhs)
{
        return (lhs / rhs.value());
}
Double_t operator / (const LauParameter& lhs, const LauParameter& rhs)
{
        return (lhs.value() / rhs.value());
}
Double_t operator += (Double_t& lhs, const LauParameter& rhs)
{
        lhs += rhs.value();
        return lhs;
}
Double_t operator -= (Double_t& lhs, const LauParameter& rhs)
{
        lhs -= rhs.value();
        return lhs;
}
Double_t operator *= (Double_t& lhs, const LauParameter& rhs)
{
        lhs *= rhs.value();
        return lhs;
}
Double_t operator /= (Double_t& lhs, const LauParameter& rhs)
{
        lhs /= rhs.value();
        return lhs;
}

// ostream operator
std::ostream& operator << (std::ostream& stream, const LauParameter& par)
{
        stream << par.value();
        return stream;
}