LauParameter.cc

Go to the documentation of this file.

/*
Copyright 2006 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 <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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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),
        secondStage_(kFALSE),
        gaussConstraint_(kFALSE),
        constraintMean_(0.0),
        constraintWidth_(0.0),
        gcc_(0.0),
        bias_(0.0),
        pull_(0.0),
        clone_(kFALSE),
        parent_(0),
        blinder_(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_),
        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_),
        blinder_((rhs.blinder_==0) ? 0 : new LauBlind(*(rhs.blinder_)))
{
}

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_;
                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_;
                delete blinder_;
                blinder_ = (rhs.blinder_==0) ? 0 : new LauBlind(*(rhs.blinder_));
        }
        return *this;
}

LauParameter::~LauParameter()
{
        delete blinder_;
}

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

void LauParameter::blindParameter(const TString& blindingString, const Double_t width)
{
        if (this->clone()) {
                parent_->blindParameter(blindingString,width);
                return;
        }

        if ( blinder_ != 0 ) {
                std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for this parameter" << std::endl;
                return;
        }

        blinder_ = new LauBlind(blindingString,width);

        for (map<LauParameter*,Double_t>::iterator iter = clones_.begin(); iter != clones_.end(); ++iter) {
                        LauParameter* clonePar = iter->first;
                        if ( clonePar->blinder_ != 0 ) {
                                std::cerr << "WARNING in LauParameter::blindParameter : blinding has already been set up for a clone of this parameter - it will be replaced!" << std::endl;
                                delete clonePar->blinder_;
                                clonePar->blinder_ = 0;
                        }
                        clonePar->blinder_ = new LauBlind(*blinder_);
        }
}

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);
        Double_t newValue = clonePar->value() * constFactor;
        clonePar->value( newValue );
        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->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);
}

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