LauCartesianGammaCPCoeffSet.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 <fstream>
#include <vector>

#include "TMath.h"
#include "TRandom.h"

#include "LauCartesianGammaCPCoeffSet.hh"
#include "LauComplex.hh"
#include "LauConstants.hh"
#include "LauParameter.hh"
#include "LauPrint.hh"

ClassImp(LauCartesianGammaCPCoeffSet)


LauCartesianGammaCPCoeffSet::LauCartesianGammaCPCoeffSet(const TString& compName, const Double_t x, const Double_t y, const Double_t xCP, const Double_t yCP, const Double_t deltaXCP, const Double_t deltaYCP,
                const Bool_t xFixed, const Bool_t yFixed, const Bool_t xCPFixed, const Bool_t yCPFixed, const Bool_t deltaXCPFixed, const Bool_t deltaYCPFixed, const Bool_t deltaXCPSecondStage, const Bool_t deltaYCPSecondStage) :
        LauAbsCoeffSet(compName),
        x_(new LauParameter("X", x, minRealImagPart_, maxRealImagPart_, xFixed)),
        y_(new LauParameter("Y", y, minRealImagPart_, maxRealImagPart_, yFixed)),
        xCP_(new LauParameter("XCP", xCP, minRealImagPart_, maxRealImagPart_, xCPFixed)),
        yCP_(new LauParameter("YCP", yCP, minRealImagPart_, maxRealImagPart_, yCPFixed)),
        deltaXCP_(new LauParameter("DeltaXCP", deltaXCP, minDelta_, maxDelta_, deltaXCPFixed)),
        deltaYCP_(new LauParameter("DeltaYCP", deltaYCP, minDelta_, maxDelta_, deltaYCPFixed)),
        nonCPPart_( x, y),
        cpPart_( 1+xCP+deltaXCP, yCP+deltaYCP),
        cpAntiPart_( 1+xCP-deltaXCP, yCP-deltaYCP),
        particleCoeff_( nonCPPart_ * cpPart_ ),
        antiparticleCoeff_( nonCPPart_ * cpAntiPart_ ),
        acp_("ACP", (antiparticleCoeff_.abs2()-particleCoeff_.abs2())/(antiparticleCoeff_.abs2()+particleCoeff_.abs2()), -1.0, 1.0, deltaXCPFixed&&deltaYCPFixed)
{
        if (deltaXCPSecondStage && !deltaXCPFixed) {
                deltaXCP_->secondStage(kTRUE);
                deltaXCP_->initValue(0.0);
        }
        if (deltaYCPSecondStage && !deltaYCPFixed) {
                deltaYCP_->secondStage(kTRUE);
                deltaYCP_->initValue(0.0);
        }
}

LauCartesianGammaCPCoeffSet::LauCartesianGammaCPCoeffSet(const LauCartesianGammaCPCoeffSet& rhs, CloneOption cloneOption, Double_t constFactor) : LauAbsCoeffSet(rhs.name()),
        x_(0),
        y_(0),
        xCP_(0),
        yCP_(0),
        deltaXCP_(0),
        deltaYCP_(0),
        nonCPPart_( rhs.nonCPPart_ ),
        cpPart_( rhs.cpPart_ ),
        cpAntiPart_( rhs.cpAntiPart_ ),
        particleCoeff_( rhs.particleCoeff_ ),
        antiparticleCoeff_( rhs.antiparticleCoeff_ ),
        acp_( rhs.acp_ )
{
        if ( cloneOption == All || cloneOption == TieRealPart ) {
                x_ = rhs.x_->createClone(constFactor);
        } else {
                x_ = new LauParameter("X", rhs.x_->value(), minRealImagPart_, maxRealImagPart_, rhs.x_->fixed());
                if ( rhs.x_->blind() ) {
                        const LauBlind* blinder = rhs.x_->blinder();
                        x_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }

        if ( cloneOption == All || cloneOption == TieImagPart ) {
                y_ = rhs.y_->createClone(constFactor);
        } else {
                y_ = new LauParameter("Y", rhs.y_->value(), minRealImagPart_, maxRealImagPart_, rhs.y_->fixed());
                if ( rhs.y_->blind() ) {
                        const LauBlind* blinder = rhs.y_->blinder();
                        y_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }

        if ( cloneOption == All || cloneOption == TieCPPars ) {
                xCP_ = rhs.xCP_->createClone(constFactor);
                yCP_ = rhs.yCP_->createClone(constFactor);
                deltaXCP_ = rhs.deltaXCP_->createClone(constFactor);
                deltaYCP_ = rhs.deltaYCP_->createClone(constFactor);
        } else {
                xCP_ = new LauParameter("XCP", rhs.xCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.xCP_->fixed());
                if ( rhs.xCP_->blind() ) {
                        const LauBlind* blinder = rhs.xCP_->blinder();
                        xCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
                yCP_ = new LauParameter("YCP", rhs.yCP_->value(), minRealImagPart_, maxRealImagPart_, rhs.yCP_->fixed());
                if ( rhs.yCP_->blind() ) {
                        const LauBlind* blinder = rhs.yCP_->blinder();
                        yCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
                deltaXCP_ = new LauParameter("DeltaXCP", rhs.deltaXCP_->value(), minDelta_, maxDelta_, rhs.deltaXCP_->fixed());
                deltaYCP_ = new LauParameter("DeltaYCP", rhs.deltaYCP_->value(), minDelta_, maxDelta_, rhs.deltaYCP_->fixed());
                if ( rhs.deltaXCP_->secondStage() && !rhs.deltaXCP_->fixed() ) {
                        deltaXCP_->secondStage(kTRUE);
                        deltaXCP_->initValue(0.0);
                }
                if ( rhs.deltaYCP_->secondStage() && !rhs.deltaYCP_->fixed() ) {
                        deltaYCP_->secondStage(kTRUE);
                        deltaYCP_->initValue(0.0);
                }
                if ( rhs.deltaXCP_->blind() ) {
                        const LauBlind* blinder = rhs.deltaXCP_->blinder();
                        deltaXCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
                if ( rhs.deltaYCP_->blind() ) {
                        const LauBlind* blinder = rhs.deltaYCP_->blinder();
                        deltaYCP_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }
}

std::vector<LauParameter*> LauCartesianGammaCPCoeffSet::getParameters()
{
        std::vector<LauParameter*> pars;
        pars.push_back(x_);
        pars.push_back(y_);
        if(!xCP_->fixed()) pars.push_back(xCP_);
        if(!yCP_->fixed()) pars.push_back(yCP_);
        if(!deltaXCP_->fixed()) pars.push_back(deltaXCP_);
        if(!deltaYCP_->fixed()) pars.push_back(deltaYCP_);
        return pars;
}

void LauCartesianGammaCPCoeffSet::printParValues() const
{
        std::cout<<"INFO in LauCartesianGammaCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
        std::cout<<"x = "<<x_->value()<<",\t";
        std::cout<<"y = "<<y_->value()<<",\t";
        std::cout<<"xCP = "<<xCP_->value()<<",\t";
        std::cout<<"yCP = "<<yCP_->value()<<",\t";
        std::cout<<"Delta xCP = "<<deltaXCP_->value()<<",\t";
        std::cout<<"Delta yCP = "<<deltaYCP_->value()<<"."<<std::endl;
}

void LauCartesianGammaCPCoeffSet::printTableHeading(std::ostream& stream) const
{
        stream<<"\\begin{tabular}{|l|c|c|c|c|c|c|}"<<std::endl;
        stream<<"\\hline"<<std::endl;
        stream<<"Component & Real Part & Imaginary Part & Real CP Part & Imaginary CP Part & $\\Delta$ Real CP Part & $\\Delta$ Imaginary CP Part \\\\"<<std::endl;
        stream<<"\\hline"<<std::endl;
}

void LauCartesianGammaCPCoeffSet::printTableRow(std::ostream& stream) const
{
        LauPrint print;
        TString resName = this->name();
        resName = resName.ReplaceAll("_", "\\_");
        stream<<resName<<"  &  $";
        print.printFormat(stream, x_->value());
        stream<<" \\pm ";
        print.printFormat(stream, x_->error());
        stream<<"$  &  $";
        print.printFormat(stream, y_->value());
        stream<<" \\pm ";
        print.printFormat(stream, y_->error());
        stream<<"$  &  $";
        print.printFormat(stream, xCP_->value());
        stream<<" \\pm ";
        print.printFormat(stream, xCP_->error());
        stream<<"$  &  $";
        print.printFormat(stream, yCP_->value());
        stream<<" \\pm ";
        print.printFormat(stream, yCP_->error());
        stream<<"$  &  $";
        print.printFormat(stream, deltaXCP_->value());
        stream<<" \\pm ";
        print.printFormat(stream, deltaXCP_->error());
        stream<<"$  &  $";
        print.printFormat(stream, deltaYCP_->value());
        stream<<" \\pm ";
        print.printFormat(stream, deltaYCP_->error());
        stream<<"$ \\\\"<<std::endl;
}

void LauCartesianGammaCPCoeffSet::randomiseInitValues()
{
        if (x_->fixed() == kFALSE) {
                // Choose a value for "X" between -3.0 and 3.0
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
                x_->initValue(value); x_->value(value);
        }
        if (y_->fixed() == kFALSE) {
                // Choose a value for "Y" between -3.0 and 3.0
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
                y_->initValue(value);  y_->value(value);
        }
        if (xCP_->fixed() == kFALSE) {
                // Choose a value for "XCP" between -3.0 and 3.0
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
                xCP_->initValue(value); xCP_->value(value);
        }
        if (yCP_->fixed() == kFALSE) {
                // Choose a value for "YCP" between -3.0 and 3.0
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*6.0 - 3.0;
                yCP_->initValue(value);  yCP_->value(value);
        }
        if (deltaXCP_->fixed() == kFALSE && deltaXCP_->secondStage() == kFALSE) {
                // Choose a value for "Delta XCP" between -0.5 and 0.5
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*1.0 - 0.5;
                deltaXCP_->initValue(value); deltaXCP_->value(value);
        }
        if (deltaYCP_->fixed() == kFALSE && deltaYCP_->secondStage() == kFALSE) {
                // Choose a value for "Delta YCP" between -0.5 and 0.5
                Double_t value = LauAbsCoeffSet::getRandomiser()->Rndm()*1.0 - 0.5;
                deltaYCP_->initValue(value); deltaYCP_->value(value);
        }
}

void LauCartesianGammaCPCoeffSet::finaliseValues()
{
        // update the pulls
        x_->updatePull();
        y_->updatePull();
        xCP_->updatePull();
        yCP_->updatePull();
        deltaXCP_->updatePull();
        deltaYCP_->updatePull();
}

const LauComplex& LauCartesianGammaCPCoeffSet::particleCoeff()
{
        nonCPPart_.setRealImagPart( x_->unblindValue(), y_->unblindValue() );
        cpPart_.setRealImagPart( 1.0+xCP_->unblindValue()+deltaXCP_->unblindValue(), yCP_->unblindValue()+deltaYCP_->unblindValue() );
        particleCoeff_ = nonCPPart_ * cpPart_;
        return particleCoeff_;
}

const LauComplex& LauCartesianGammaCPCoeffSet::antiparticleCoeff()
{
        nonCPPart_.setRealImagPart( x_->unblindValue(), y_->unblindValue() );
        cpAntiPart_.setRealImagPart( 1.0+xCP_->unblindValue()-deltaXCP_->unblindValue(), yCP_->unblindValue()-deltaYCP_->unblindValue() );
        antiparticleCoeff_ = nonCPPart_ * cpAntiPart_;
        return antiparticleCoeff_;
}

void LauCartesianGammaCPCoeffSet::setCoeffValues( const LauComplex&, const LauComplex&, Bool_t )
{
        std::cerr << "ERROR in LauCartesianGammaCPCoeffSet::setCoeffValues : Method not supported by this class - too many parameters" << std::endl;
}

LauParameter LauCartesianGammaCPCoeffSet::acp()
{
        // set the name
        TString parName(this->baseName()); parName += "_ACP";
        acp_.name(parName);

        // work out the ACP value
        const LauComplex nonCPPart( x_->value(), y_->value() );
        const LauComplex cpPart( 1.0+xCP_->value()+deltaXCP_->value(), yCP_->value()+deltaYCP_->value() );
        const LauComplex cpAntiPart( 1.0+xCP_->value()-deltaXCP_->value(), yCP_->value()-deltaYCP_->value() );
        const LauComplex partCoeff = nonCPPart * cpPart;
        const LauComplex antiCoeff = nonCPPart * cpAntiPart;

        const Double_t numer = antiCoeff.abs2() - partCoeff.abs2();
        const Double_t denom = antiCoeff.abs2() + partCoeff.abs2();
        const Double_t value = numer/denom;

        // is it fixed?
        const Bool_t fixed = deltaXCP_->fixed() && deltaYCP_->fixed();
        acp_.fixed(fixed);

        // we can't work out the error without the covariance matrix
        const Double_t error(0.0);

        // set the value and error
        acp_.valueAndErrors(value,error);

        return acp_;
}

LauAbsCoeffSet* LauCartesianGammaCPCoeffSet::createClone(const TString& newName, CloneOption cloneOption, Double_t constFactor)
{
        LauAbsCoeffSet* clone(0);
        if ( cloneOption == All || cloneOption == TieRealPart || cloneOption == TieImagPart || cloneOption == TieCPPars ) {
                clone = new LauCartesianGammaCPCoeffSet( *this, cloneOption, constFactor );
                clone->name( newName );
        } else {
                std::cerr << "ERROR in LauCartesianGammaCPCoeffSet::createClone : Invalid clone option" << std::endl;
        }
        return clone;
}