LauBelleCPCoeffSet.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 "LauComplex.hh"
#include "LauConstants.hh"
#include "LauBelleCPCoeffSet.hh"
#include "LauParameter.hh"
#include "LauPrint.hh"
#include "LauRandom.hh"

ClassImp(LauBelleCPCoeffSet)


LauBelleCPCoeffSet::LauBelleCPCoeffSet(const TString& compName, Double_t a, Double_t delta, Double_t b, Double_t phi,
                Bool_t aFixed, Bool_t deltaFixed, Bool_t bFixed, Bool_t phiFixed, Bool_t bSecondStage, Bool_t phiSecondStage) :
        LauAbsCoeffSet(compName),
        a_(new LauParameter("A", a, minMagnitude_, maxMagnitude_, aFixed)),
        b_(new LauParameter("B", b, minMagnitude_, maxMagnitude_, bFixed)),
        delta_(new LauParameter("Delta", delta, minPhase_, maxPhase_, deltaFixed)),
        phi_(new LauParameter("Phi", phi, minPhase_, maxPhase_, phiFixed)),
        particleCoeff_(0.0,0.0),
        antiparticleCoeff_(0.0,0.0),
        acp_("ACP", (-2.0*b*TMath::Cos(phi))/(1.0+b*b), -1.0, 1.0, bFixed&&phiFixed)
{
        if (bSecondStage && !bFixed) {
                b_->secondStage(kTRUE);
                b_->initValue(0.0);
        }
        if (phiSecondStage && !phiFixed) {
                phi_->secondStage(kTRUE);
                phi_->initValue(0.0);
        }
}

LauBelleCPCoeffSet::LauBelleCPCoeffSet(const LauBelleCPCoeffSet& rhs, CloneOption cloneOption, Double_t constFactor) : LauAbsCoeffSet(rhs.name()),
        a_(0),
        b_(0),
        delta_(0),
        phi_(0),
        particleCoeff_( rhs.particleCoeff_ ),
        antiparticleCoeff_( rhs.antiparticleCoeff_ ),
        acp_( rhs.acp_ )
{
        if ( cloneOption == All || cloneOption == TieMagnitude ) {
                a_ = rhs.a_->createClone(constFactor);
        } else {
                a_ = new LauParameter("A", rhs.a_->value(), minMagnitude_, maxMagnitude_, rhs.a_->fixed());
                if ( rhs.a_->blind() ) {
                        const LauBlind* blinder = rhs.a_->blinder();
                        a_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }

        if ( cloneOption == All || cloneOption == TieCPPars ) {
                b_ = rhs.b_->createClone(constFactor);
        } else {
                b_ = new LauParameter("B", rhs.b_->value(), minMagnitude_, maxMagnitude_, rhs.b_->fixed());
                if ( rhs.b_->blind() ) {
                        const LauBlind* blinder = rhs.b_->blinder();
                        b_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }

        if ( cloneOption == All || cloneOption == TiePhase ) {
                delta_ = rhs.delta_->createClone(constFactor);
        } else {
                delta_ = new LauParameter("Delta", rhs.delta_->value(), minPhase_, maxPhase_, rhs.delta_->fixed());
                if ( rhs.delta_->blind() ) {
                        const LauBlind* blinder = rhs.delta_->blinder();
                        delta_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }

        if ( cloneOption == All || cloneOption == TieCPPars ) {
                phi_ = rhs.phi_->createClone(constFactor);
        } else {
                phi_ = new LauParameter("Phi", rhs.phi_->value(), minPhase_, maxPhase_, rhs.phi_->fixed());
                if ( rhs.phi_->blind() ) {
                        const LauBlind* blinder = rhs.phi_->blinder();
                        phi_->blindParameter( blinder->blindingString(), blinder->blindingWidth() );
                }
        }
}

std::vector<LauParameter*> LauBelleCPCoeffSet::getParameters()
{
        std::vector<LauParameter*> pars;
        pars.push_back(a_);
        pars.push_back(b_);
        pars.push_back(delta_);
        pars.push_back(phi_);
        return pars;
}

void LauBelleCPCoeffSet::printParValues() const
{
        std::cout<<"INFO in LauBelleCPCoeffSet::printParValues : Component \""<<this->name()<<"\" has ";
        std::cout<<"a-magnitude = "<<a_->value()<<",\t";
        std::cout<<"delta = "<<delta_->value()<<",\t";
        std::cout<<"b-magnitude = "<<b_->value()<<",\t";
        std::cout<<"phi = "<<phi_->value()<<"."<<std::endl;
}

void LauBelleCPCoeffSet::printTableHeading(std::ostream& stream) const
{
        stream<<"\\begin{tabular}{|l|c|c|c|c|}"<<std::endl;
        stream<<"\\hline"<<std::endl;
        stream<<"Component & a-Magnitude & delta & b-Magnitude & phi \\\\"<<std::endl;
        stream<<"\\hline"<<std::endl;
}

void LauBelleCPCoeffSet::printTableRow(std::ostream& stream) const
{
        LauPrint print;
        TString resName = this->name();
        resName = resName.ReplaceAll("_", "\\_");
        stream<<resName<<"  &  $";
        print.printFormat(stream, a_->value());
        stream<<" \\pm ";
        print.printFormat(stream, a_->error());
        stream<<"$  &  $";
        print.printFormat(stream, delta_->value());
        stream<<" \\pm ";
        print.printFormat(stream, delta_->error());
        stream<<"$  &  $";
        print.printFormat(stream, b_->value());
        stream<<" \\pm ";
        print.printFormat(stream, b_->error());
        stream<<"$  &  $";
        print.printFormat(stream, phi_->value());
        stream<<" \\pm ";
        print.printFormat(stream, phi_->error());
        stream<<"$ \\\\"<<std::endl;
}

void LauBelleCPCoeffSet::randomiseInitValues()
{
        if (a_->fixed() == kFALSE) {
                // Choose an a-magnitude between 0.0 and 2.0
                Double_t mag = LauRandom::zeroSeedRandom()->Rndm()*2.0;
                a_->initValue(mag); a_->value(mag);
        }
        if (b_->fixed() == kFALSE && b_->secondStage() == kFALSE) {
                // Choose a b-magnitude between 0.0 and 0.1
                Double_t mag = LauRandom::zeroSeedRandom()->Rndm()*0.1;
                b_->initValue(mag);  b_->value(mag);
        }
        if (delta_->fixed() == kFALSE) {
                // Choose a phase between +- pi
                Double_t phase = LauRandom::zeroSeedRandom()->Rndm()*LauConstants::twoPi - LauConstants::pi;
                delta_->initValue(phase); delta_->value(phase);
        }
        if (phi_->fixed() == kFALSE && phi_->secondStage() == kFALSE) {
                // Choose a phase between +- pi
                Double_t phase = LauRandom::zeroSeedRandom()->Rndm()*LauConstants::twoPi - LauConstants::pi;
                phi_->initValue(phase); phi_->value(phase);
        }
}

void LauBelleCPCoeffSet::finaliseValues()
{
        // retrieve the current values from the parameters
        Double_t aVal     = a_->value();
        Double_t bVal     = b_->value();
        Double_t deltaVal = delta_->value();
        Double_t phiVal   = phi_->value();
        Double_t genDelta = delta_->genValue();
        Double_t genPhi   = phi_->genValue();

        // Check whether we have a negative "a" magnitude.
        // If so make it positive and add pi to the "delta" phase.
        if (aVal < 0.0) {
                aVal *= -1.0;
                deltaVal += LauConstants::pi;
        }

        // Check whether we have a negative "b" magnitude.
        // If so make it positive and add pi to the "phi" phase.
        if (bVal < 0.0) {
                bVal *= -1.0;
                phiVal += LauConstants::pi;
        }

        // Check now whether the phases lies in the right range (-pi to pi).
        Bool_t deltaWithinRange(kFALSE);
        Bool_t phiWithinRange(kFALSE);
        while (deltaWithinRange == kFALSE && phiWithinRange == kFALSE) {
                if (deltaVal > -LauConstants::pi && deltaVal < LauConstants::pi) {
                        deltaWithinRange = kTRUE;
                } else {
                        // Not within the specified range
                        if (deltaVal > LauConstants::pi) {
                                deltaVal -= LauConstants::twoPi;
                        } else if (deltaVal < -LauConstants::pi) {
                                deltaVal += LauConstants::twoPi;
                        }
                }

                if (phiVal > -LauConstants::pi && phiVal < LauConstants::pi) {
                        phiWithinRange = kTRUE;
                } else {
                        // Not within the specified range
                        if (phiVal > LauConstants::pi) {
                                phiVal -= LauConstants::twoPi;
                        } else if (phiVal < -LauConstants::pi) {
                                phiVal += LauConstants::twoPi;
                        }
                }
        }

        // A further problem can occur when the generated phase is close to -pi or pi.
        // The phase can wrap over to the other end of the scale -
        // this leads to artificially large pulls so we wrap it back.
        Double_t diff = deltaVal - genDelta;
        if (diff > LauConstants::pi) {
                deltaVal -= LauConstants::twoPi;
        } else if (diff < -LauConstants::pi) {
                deltaVal += LauConstants::twoPi;
        }

        diff = phiVal - genPhi;
        if (diff > LauConstants::pi) {
                phiVal -= LauConstants::twoPi;
        } else if (diff < -LauConstants::pi) {
                phiVal += LauConstants::twoPi;
        }

        // finally store the new values in the parameters
        // and update the pulls
        a_->value(aVal);          a_->updatePull();
        b_->value(bVal);          b_->updatePull();
        delta_->value(deltaVal);  delta_->updatePull();
        phi_->value(phiVal);      phi_->updatePull();
}

const LauComplex& LauBelleCPCoeffSet::particleCoeff()
{
        LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
        LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
        particleCoeff_.setRealImagPart(1.0,0.0);
        particleCoeff_ += bTerm;
        particleCoeff_ *= aTerm;
        return particleCoeff_;
}

const LauComplex& LauBelleCPCoeffSet::antiparticleCoeff()
{
        LauComplex aTerm(a_->unblindValue()*TMath::Cos(delta_->unblindValue()), a_->unblindValue()*TMath::Sin(delta_->unblindValue()));
        LauComplex bTerm(b_->unblindValue()*TMath::Cos(phi_->unblindValue()), b_->unblindValue()*TMath::Sin(phi_->unblindValue()));
        antiparticleCoeff_.setRealImagPart(1.0,0.0);
        antiparticleCoeff_ -= bTerm;
        antiparticleCoeff_ *= aTerm;
        return antiparticleCoeff_;
}

void LauBelleCPCoeffSet::setCoeffValues( const LauComplex& coeff, const LauComplex& coeffBar, Bool_t init )
{
        LauComplex sum = coeff + coeffBar;
        LauComplex diff = coeff - coeffBar;
        LauComplex ratio = diff / sum;

        Double_t aVal( 0.5 * sum.abs() );
        Double_t deltaVal( sum.arg() );
        Double_t bVal( ratio.abs() );
        Double_t phiVal( ratio.arg() );

        a_->value( aVal );
        delta_->value( deltaVal );
        b_->value( bVal );
        phi_->value( phiVal );

        if ( init ) {
                a_->genValue( aVal );
                delta_->genValue( deltaVal );
                b_->genValue( bVal );
                phi_->genValue( phiVal );

                a_->initValue( aVal );
                delta_->initValue( deltaVal );
                b_->initValue( bVal );
                phi_->initValue( phiVal );
        }
}

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

        // work out the ACP value
        Double_t value = (-2.0*b_->value()*TMath::Cos(phi_->value()))/(1.0+b_->value()*b_->value());

        // is it fixed?
        Bool_t fixed = b_->fixed() && phi_->fixed();
        acp_.fixed(fixed);

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

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

        return acp_;
}

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