chisquare.cpp File Reference

#include <goofit/Application.h>
#include <goofit/BinnedDataSet.h>
#include <goofit/FitControl.h>
#include <goofit/FitManager.h>
#include <goofit/UnbinnedDataSet.h>
#include <goofit/Variable.h>
#include <CLI/Timer.hpp>
#include <TCanvas.h>
#include <TH1F.h>
#include <TLatex.h>
#include <TMinuit.h>
#include <TRandom.h>
#include <goofit/PDFs/basic/PolynomialPdf.h>
#include <iostream>
#include <string>
#include <tuple>
#include <vector>

Go to the source code of this file.

Functions
double	integralExpCon (double lo, double hi)
double	integralExpLin (double lo, double hi)
double	integralExpSqu (double lo, double hi)
void	generateEvents (Observable decayTime, vector< int > &rsEvtVec, vector< int > &wsEvtVec, double conCoef, double linCoef, double squCoef, int eventsToGenerate)
std::tuple< int, std::string >	fitRatio (Observable decayTime, vector< Variable > weights, vector< int > &rsEvts, vector< int > &wsEvts, std::string plotName="")
void	cpvFitFcn (int &npar, double *gin, double &fun, double *fp, int iflag)
void	fitRatioCPU (Observable decayTime, vector< int > &rsEvts, vector< int > &wsEvts)
int	main (int argc, char **argv)

Variables
TCanvas	foo
vector< double >	ratios
vector< double >	errors
Observable	decayTime {"decayTime", 0, 10}

Function Documentation

cpvFitFcn()

void cpvFitFcn	(	int &	npar,
		double *	gin,
		double &	fun,
		double *	fp,
		int	iflag
	)

Definition at line 160 of file chisquare.cpp.

References decayTime, errors, GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::Indexable::getUpperLimit(), and ratios.

Referenced by fitRatioCPU().

                                                                           {
    double conCoef = fp[0];
    double linCoef = fp[1];
    double squCoef = fp[2];

    double chisq = 0;
    double step  = (decayTime.getUpperLimit() - decayTime.getLowerLimit()) / decayTime.getNumBins();

    for(unsigned int i = 0; i < ratios.size(); ++i) {
        double currDTime = decayTime.getLowerLimit() + (i + 0.5) * step;
        double pdfval    = conCoef + linCoef * currDTime + squCoef * currDTime * currDTime;
        chisq += pow((pdfval - ratios[i]) / errors[i], 2);
    }

    fun = chisq;
}

fitRatio()

std::tuple<int, std::string> fitRatio	(	Observable	decayTime,
		vector< Variable >	weights,
		vector< int > &	rsEvts,
		vector< int > &	wsEvts,
		std::string	plotName = ""
	)

Definition at line 76 of file chisquare.cpp.

References GooFit::GooPdf::evaluateAtPoints(), GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::Indexable::getUpperLimit(), GooFit::PdfBase::setData(), and GooFit::GooPdf::setFitControl().

Referenced by main().

                                                               {
    TH1D ratioHist("ratioHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());

    BinnedDataSet ratioData(decayTime);

    for(unsigned int i = 0; i < wsEvts.size(); ++i) {
        double ratio = wsEvts[i];

        if(0 == rsEvts[i])
            rsEvts[i] = 1; // Cheating to avoid div by zero.

        ratio /= rsEvts[i];

        if(0 == wsEvts[i])
            wsEvts[i] = 1; // Avoid zero errors

        double error = wsEvts[i] / pow(rsEvts[i], 2);
        error += pow(wsEvts[i], 2) / pow(rsEvts[i], 3);
        error = sqrt(error);

        ratioData.setBinContent(i, ratio);
        ratioData.setBinError(i, error);
        ratioHist.SetBinContent(i + 1, ratio);
        ratioHist.SetBinError(i + 1, error);
    }

    PolynomialPdf poly("poly", decayTime, weights);
    poly.setFitControl(std::make_shared<BinnedErrorFit>());
    poly.setData(&ratioData);
    FitManager fitter{&poly};

    CLI::Timer timer_cpu{"GPU"};
    fitter.fit();
    std::string timer_str = timer_cpu.to_string();

    if(!plotName.empty()) {
        vector<fptype> values = poly.evaluateAtPoints(decayTime);
        TH1D pdfHist("pdfHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());

        for(int i = 0; i < values.size(); ++i) {
            pdfHist.SetBinContent(i + 1, values[i]);
        }

        ratioHist.SetMarkerStyle(8);
        ratioHist.SetMarkerSize(0.5);
        ratioHist.SetStats(false);
        ratioHist.Draw("p");

        TString str1 = fmt::format(
            "Constant [10^{{-2}}] : {:.3} #pm {:.3}", weights[0].getValue() * 1e2, weights[0].getError() * 1e2);
        TLatex res1(0.14, 0.83, str1);
        res1.SetNDC(true);

        TString str2 = fmt::format(
            "Linear [10^{{-4}}] : {:.3} #pm {:.3}", weights[1].getValue() * 1e4, weights[1].getError() * 1e4);
        TLatex res2(0.14, 0.73, str2);
        res2.SetNDC(true);

        TString str3 = fmt::format(
            "Quadratic [10^{{-6}}]: {:.3} #pm {:.3}", weights[2].getValue() * 1e6, weights[2].getError() * 1e6);
        TLatex res3(0.14, 0.63, str3);
        res3.SetNDC(true);

        res1.Draw();
        res2.Draw();
        res3.Draw();

        pdfHist.SetLineColor(kBlue);
        pdfHist.SetLineWidth(3);
        pdfHist.SetStats(false);
        pdfHist.Draw("lsame");
        foo.SaveAs(plotName.c_str());
    }

    std::cout << "Polynomial function: " << poly.getCoefficient(2) << " * t^2 + " << poly.getCoefficient(1) << " * t + "
              << poly.getCoefficient(0) << std::endl;

    return make_tuple(int(fitter), timer_str);
}

fitRatioCPU()

void fitRatioCPU	(	Observable	decayTime,
		vector< int > &	rsEvts,
		vector< int > &	wsEvts
	)

Definition at line 177 of file chisquare.cpp.

References cpvFitFcn(), errors, GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::Indexable::getUpperLimit(), and ratios.

Referenced by main().

                                                                                 {
    TH1D *ratioHist
        = new TH1D("ratioHist", "", decayTime.getNumBins(), decayTime.getLowerLimit(), decayTime.getUpperLimit());

    ratios.resize(wsEvts.size());
    errors.resize(wsEvts.size());

    for(unsigned int i = 0; i < wsEvts.size(); ++i) {
        if(0 == rsEvts[i])
            rsEvts[i] = 1; // Cheating to avoid div by zero.

        fptype ratio = wsEvts[i] / rsEvts[i];

        if(0 == wsEvts[i])
            wsEvts[i] = 1; // Avoid zero errors

        double error = wsEvts[i] / pow(rsEvts[i], 2);
        error += pow(wsEvts[i], 2) / pow(rsEvts[i], 3);
        error = sqrt(error);

        ratios[i] = ratio;
        errors[i] = error;
        ratioHist->SetBinContent(i + 1, ratio);
        ratioHist->SetBinError(i + 1, error);
    }

    TMinuit *minuit = new TMinuit(3);
    minuit->DefineParameter(0, "constaCoef", 0.03, 0.01, -1, 1);
    minuit->DefineParameter(1, "linearCoef", 0, 0.01, -1, 1);
    minuit->DefineParameter(2, "secondCoef", 0, 0.01, -1, 1);
    minuit->SetFCN(cpvFitFcn);

    minuit->Migrad();
}

generateEvents()

void generateEvents	(	Observable	decayTime,
		vector< int > &	rsEvtVec,
		vector< int > &	wsEvtVec,
		double	conCoef,
		double	linCoef,
		double	squCoef,
		int	eventsToGenerate
	)

Definition at line 41 of file chisquare.cpp.

References GooFit::Indexable::getLowerLimit(), GooFit::Observable::getNumBins(), GooFit::Indexable::getUpperLimit(), integralExpCon(), integralExpLin(), and integralExpSqu().

Referenced by main().

                                          {
    static TRandom donram(24);
    double totalRSintegral = integralExpCon(0, 100);
    double step            = (decayTime.getUpperLimit() - decayTime.getLowerLimit()) / decayTime.getNumBins();

    for(int i = 0; i < decayTime.getNumBins(); ++i) {
        double binStart = i * step;
        binStart += decayTime.getLowerLimit();
        double binFinal = binStart + step;

        double rsIntegral = integralExpCon(binStart, binFinal);
        double wsIntegral = conCoef * integralExpCon(binStart, binFinal);
        wsIntegral += linCoef * integralExpLin(binStart, binFinal);
        wsIntegral += squCoef * integralExpSqu(binStart, binFinal);

        double expectedRSevts = eventsToGenerate * rsIntegral / totalRSintegral;
        double expectedWSevts = eventsToGenerate * wsIntegral / totalRSintegral;

        int rsEvts  = donram.Poisson(expectedRSevts);
        int wsEvts  = donram.Poisson(expectedWSevts);
        rsEvtVec[i] = rsEvts;
        wsEvtVec[i] = wsEvts;

        if(0 == (i % 10))
            std::cout << "Events in bin " << i << " : " << rsEvts << " (" << expectedRSevts << ") " << wsEvts << " ("
                      << expectedWSevts << ")\n";
    }
}

integralExpCon()

double integralExpCon	(	double	lo,
		double	hi
	)

Definition at line 33 of file chisquare.cpp.

Referenced by generateEvents().

{ return (exp(-lo) - exp(-hi)); }

integralExpLin()

double integralExpLin	(	double	lo,
		double	hi
	)

Definition at line 35 of file chisquare.cpp.

Referenced by generateEvents().

{ return ((lo + 1) * exp(-lo) - (hi + 1) * exp(-hi)); }

integralExpSqu()

double integralExpSqu	(	double	lo,
		double	hi
	)

Definition at line 37 of file chisquare.cpp.

Referenced by generateEvents().

                                            {
    return ((lo * lo + 2 * lo + 2) * exp(-lo) - (hi * hi + 2 * hi + 2) * exp(-hi));
}

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 212 of file chisquare.cpp.

References decayTime, fitRatio(), fitRatioCPU(), generateEvents(), GooFit::Observable::getNumBins(), GOOFIT_PARSE, GooFit::Observable::setNumBins(), GooFit::Indexable::setValue(), and weights.

                                {
    GooFit::Application app("Chi-square example", argc, argv);

    int numbins = 100;
    app.add_option("-n,--numbins", numbins, "Number of bins", true);

    int eventsToGenerate = 10000000;
    app.add_option("-e,--events", eventsToGenerate, "Events to generate", true);

    GOOFIT_PARSE(app);

    // Time is in units of lifetime

    decayTime.setValue(100.);
    decayTime.setNumBins(numbins);

    double rSubD = 0.03;
    double rBarD = 0.03;
    double delta = 0;
    double wpPhi = 0;
    double x_mix = 0.0016;
    double y_mix = 0.0055;
    double magPQ = 1.0;
    double magQP = 1.0 / magPQ;

    vector<int> dZeroEvtsWS(decayTime.getNumBins());
    vector<int> dZeroEvtsRS(decayTime.getNumBins());
    vector<int> d0barEvtsWS(decayTime.getNumBins());
    vector<int> d0barEvtsRS(decayTime.getNumBins());

    double dZeroLinearCoef = magPQ * sqrt(rSubD) * (y_mix * cos(delta + wpPhi) - x_mix * sin(delta + wpPhi));
    double d0barLinearCoef = magQP * sqrt(rBarD) * (y_mix * cos(delta - wpPhi) - x_mix * sin(delta - wpPhi));

    double dZeroSecondCoef = 0.25 * magPQ * magPQ * (x_mix * x_mix + y_mix * y_mix);
    double d0barSecondCoef = 0.25 * magQP * magQP * (x_mix * x_mix + y_mix * y_mix);

    generateEvents(decayTime, dZeroEvtsRS, dZeroEvtsWS, rSubD, dZeroLinearCoef, dZeroSecondCoef, eventsToGenerate);
    generateEvents(decayTime, d0barEvtsRS, d0barEvtsWS, rBarD, d0barLinearCoef, d0barSecondCoef, eventsToGenerate);

    Variable constaCoef("constaCoef", 0.03, 0.01, -1, 1);
    Variable linearCoef("linearCoef", 0, 0.01, -1, 1);
    Variable secondCoef("secondCoef", 0, 0.01, -1, 1);

    vector<Variable> weights = {constaCoef, linearCoef, secondCoef};

    int retval1, retval2;
    std::string fit1, fit2;
    std::tie(retval1, fit1)
        = fitRatio(decayTime, weights, dZeroEvtsRS, dZeroEvtsWS, "chisquare_dzeroEvtRatio_goo_cpp.png");
    std::tie(retval2, fit2)
        = fitRatio(decayTime, weights, d0barEvtsRS, d0barEvtsWS, "chisquare_dzbarEvtRatio_goo_cpp.png");

    CLI::Timer timer_cpu{"Total CPU (2x fits)"};
    fitRatioCPU(decayTime, dZeroEvtsRS, dZeroEvtsWS);
    fitRatioCPU(decayTime, d0barEvtsRS, d0barEvtsWS);
    std::string cpu_string = timer_cpu.to_string();

    std::cout << fit1 << "\n" << fit2 << "\n" << cpu_string << std::endl;

    fmt::print("Exit codes (should be 0): {} and {}\n", retval1, retval2);

    return retval1 + retval2;
}

Variable Documentation

decayTime

Observable decayTime {"decayTime", 0, 10}

Definition at line 31 of file chisquare.cpp.

Referenced by cpvFitFcn(), and main().

errors

vector<double> errors

Definition at line 28 of file chisquare.cpp.

Referenced by cpvFitFcn(), and fitRatioCPU().

foo

TCanvas foo

Definition at line 22 of file chisquare.cpp.

Referenced by getToyData(), main(), and runToyFit().

ratios

vector<double> ratios

Definition at line 27 of file chisquare.cpp.

Referenced by cpvFitFcn(), and fitRatioCPU().