//Renata Kopecna


#include "tests.hh"

#include <spdlog.h>
#include <integrals.hh>

//TODO: go through what is needed
#include <complex>
#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include <TMath.h>
#include <TLegend.h>
#include <TH2D.h>
#include <TRandom3.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TFile.h>
#include <TTree.h>
#include <Math/SpecFuncMathCore.h>
#include <Math/SpecFuncMathMore.h>
#include <Math/Vector3D.h>
#include <Math/Boost.h>
#include <TMatrixDSym.h>
#include <TVectorD.h>
#include <TDecompChol.h>
#include <TDecompBK.h>
#include <TMatrixD.h>
#include <Math/GSLIntegrator.h>
#include <funcs.hh>
#include <event.hh>

namespace fcnc {

void test_integrals(){
    spdlog::debug("Starting Test");
    //sin test
    unsigned int nbins=100;
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_sin_x_2(x, 0);
        double b = 0.5*(x-0.5*sin(2.0*x));
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_sin_x_2(x, 1);
        double b = -(2*x*sin(2*x)+cos(2*x)-2*x*x)/8;
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_sin_x_2(x, 2);
        double b = -((6*x*x-3)*sin(2*x)+6*x*cos(2*x)-4*x*x*x)/24;
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_sin_x_2(x, 3);
        double b = -((4*x*x*x-6*x)*sin(2*x)+(6*x*x-3)*cos(2*x)-2*x*x*x*x)/16;
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_sin_x_2(x, 1);
        double b = -(2*x*sin(2*x)+cos(2*x)-2*x*x)/8;
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    //cos test
    for (unsigned int i = 0; i < nbins; i++)
    {
        double x = -1.0 + 2.0*i/double(nbins);
        double a = integral_x_to_n_times_cos_x_2(x, 0);
        double b = (sin(2.0*x)*0.5+x)*0.5;
        if (fabs(a-b) > 1.0e-8)
        {
            std::cout << a;
            std::cout << " != ";
            std::cout << b;
            std::cout << std::endl;
        }
    }
    std::vector<double> res;
    chebyshev(0.43, 13, res);
    for (unsigned int i = 0; i < 14; i++)
    {
        if (fabs(res.at(i) - chebyshev(0.43,i)) > 1.0e-8){
            spdlog::debug("CHEBYSHEV: {0:f} != {1:f}", res.at(i), chebyshev(0.43,i));
        }
    }
    legendre(0.43, 13, res);
    for (unsigned int i = 0; i < 14; i++)
    {
        if (fabs(res.at(i) - ROOT::Math::legendre(i, 0.43)) > 1.0e-8){
            spdlog::debug("LEGENDREPOLY {0:d}: {1:f} != {2:f}", i, res.at(i), ROOT::Math::legendre(i, 0.43));
        }
    }

    spdlog::debug("Finished Test");
}

void test_legendre()
{
    spdlog::debug("Starting Test");
    TRandom3 rnd;
    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200);
    canvas->cd();
    unsigned int nbins = 200;
    TH1D * data = new TH1D("data", "data", nbins, -1.0, 1.0);
    for (unsigned int i =0; i<200000; i++){
        data->Fill(rnd.Gaus(0.0,0.5));
        bool finished = false;
        double x;
        while (!finished){
            x = rnd.Rndm();
            if (rnd.Rndm() < x){
                x = x*2.0-1.0;
                finished = true;
            }
        }
        data->Fill(x);
    }
    //data->FillRandom("gaus", 200000);
    data->Draw("e");

    unsigned int n = 20;

    std::vector<double> l;
    double dv = 2.0/nbins;
    std::vector<double> coeffs(n+1, 0.0);
    for (unsigned int j=0; j<nbins; j++) {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        //orthonormal_legendre(x, n, l);
        legendre(x, n, l);
        double y = data->GetBinContent(j+1);
        for (unsigned int i = 0; i<n+1; i++){
            coeffs.at(i) += y*l.at(i) * dv * (0.5*(2.0*i+1.0));
        }
    }
    for (unsigned int i = 0; i<n+1; i++){
        spdlog::debug("coefficient {0:d} = {1:f}", i, coeffs.at(i));
    }
    TH1D * interpol = new TH1D("inter", "inter", nbins, -1.0, 1.0);
    for (unsigned int j=0; j<nbins; j++){
        double x = -1.0 + (j+0.5)*2.0/nbins;
        //orthonormal_legendre(x, n, l);
        legendre(x, n, l);
        double y = 0.0;
        for (unsigned int i = 0; i<n+1; i++){
            y += coeffs.at(i) * l.at(i);
        }
        interpol->SetBinContent(j+1, y);
    }
    interpol->SetLineColor(2);
    interpol->Draw("same l");

    TH1D * pol1d = new TH1D("pol1d", "pol1d", nbins, -1.0, 1.0);
    std::vector<double> poly_coeffs;
    for (unsigned int i=0; i<n+1; i++){
        spdlog::debug("Poly Coefficient {0:d} = {1:f}", i, poly_coeffs.at(i));
    }
    for (unsigned int j=0; j<nbins; j++)    {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        double y = 0.0;
        for (unsigned int i = 0; i<n+1; i++){
            y += poly_coeffs.at(i) * pow(x, int(i));
        }
        pol1d->SetBinContent(j+1, y);
    }
    pol1d->SetLineColor(3);
    pol1d->SetLineStyle(2);
    pol1d->Draw("same l");

    canvas->Update();
    canvas->Print("legendre.eps", "eps");
    delete data;
    delete interpol;
    delete pol1d;
    delete canvas;
}

void test_legendre_2pi()
{
    TRandom3 rnd;
    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200) ;
    canvas->cd();
    unsigned int nbins = 200;
    TH1D * data = new TH1D("data", "data", nbins, -TMath::Pi(), TMath::Pi());
    for (unsigned int i =0; i<200000; i++)
    {
        data->Fill(rnd.Gaus(0.0,0.5));
        bool finished = false;
        double x;
        while (!finished)
        {
            x = rnd.Rndm();
            if (rnd.Rndm() < x)
            {
                x = x*2.0*TMath::Pi()-TMath::Pi();
                finished = true;
            }
        }
        data->Fill(x);
    }
    //data->FillRandom("gaus", 200000);
    data->Draw("e");

    unsigned int n = 20;

    std::vector<double> l;
    double dv = 2.0/nbins;
    std::vector<double> coeffs(n+1, 0.0);
    for (unsigned int j=0; j<nbins; j++)
    {
        double x = -1.0 + (j+0.5)*2.0/nbins; //this transforms correctly to -1...+1
        //orthonormal_legendre(x, n, l);
        legendre(x, n, l);
        double y = data->GetBinContent(j+1);
        for (unsigned int i = 0; i<n+1; i++)
            coeffs.at(i) += y*l.at(i) * dv * (0.5*(2.0*i+1.0));
    }
    for (unsigned int i = 0; i<n+1; i++){
        spdlog::debug("coefficient {0:d} = {1:f}", i, coeffs.at(i));
    }
    TH1D * interpol = new TH1D("inter", "inter", nbins, -TMath::Pi(), +TMath::Pi());
    for (unsigned int j=0; j<nbins; j++)
    {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        //double x = (-1.0 + (j+0.5)*2.0/nbins) / TMath::Pi();

        //double x = -TMath::Pi() + (j+0.5)*2.0*TMath::Pi()/nbins;
        //orthonormal_legendre(x, n, l);
        legendre(x, n, l);
        double y = 0.0;
        for (unsigned int i = 0; i<n+1; i++){
            y += coeffs.at(i) * l.at(i);
        }
        interpol->SetBinContent(j+1, y);
    }
    interpol->SetLineColor(2);
    interpol->Draw("same l");

    TH1D * pol1d = new TH1D("pol1d", "pol1d", nbins, -TMath::Pi(), +TMath::Pi());
    std::vector<double> poly_coeffs;
    for (unsigned int i=0; i<n+1; i++){
        spdlog::debug("Poly Coefficient {0:d} = {1:f}", i,  poly_coeffs.at(i));
    }
    for (unsigned int j=0; j<nbins; j++)
    {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        //double x = (-1.0 + (j+0.5)*2.0/nbins)/TMath::Pi();
        //double x = xc/TMath::Pi();
        double y = 0.0;
        for (unsigned int i = 0; i<n+1; i++){
            y += poly_coeffs.at(i) * pow(x, int(i));
        }
        pol1d->SetBinContent(j+1, y);
    }
    pol1d->SetLineColor(3);
    pol1d->SetLineStyle(2);
    pol1d->Draw("same l");

    canvas->Update();
    canvas->Print("legendre2pi.eps", "eps");
    delete data;
    delete interpol;
    delete pol1d;
    delete canvas;
}

void test_chebyshev()
{
    TRandom3 rnd;
    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200) ;
    canvas->cd();
    unsigned int nbins = 200;
    TH1D * data = new TH1D("data", "data", nbins, -1.0, 1.0);
    for (unsigned int i =0; i<200000; i++)
    {
        data->Fill(rnd.Gaus(0.0,0.5));
        bool finished = false;
        double x;
        while (!finished)
        {
            x = rnd.Rndm();
            if (rnd.Rndm() < x)
            {
                x = x*2.0-1.0;
                finished = true;
            }
        }
        data->Fill(x);
    }
    //data->FillRandom("gaus", 200000);
    data->Draw("e");

    unsigned int n = 20;

    std::vector<double> l;
    double dv = 2.0/nbins;
    std::vector<double> coeffs(n+1, 0.0);
    for (unsigned int j=0; j<nbins; j++)
    {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        chebyshev(x, n, l);
        double y = data->GetBinContent(j+1);
        for (unsigned int i = 0; i<n+1; i++){
            coeffs.at(i) += y*l.at(i) * dv / sqrt(1-x*x) / (i==0 ? TMath::Pi() : 0.5*TMath::Pi() );
        }
    }
    for (unsigned int i = 0; i<n+1; i++){
        spdlog::debug("coefficient {0:d} = {1:f}", i,  coeffs.at(i));
    }
    TH1D * interpol = new TH1D("inter", "inter", nbins, -1.0, 1.0);
    for (unsigned int j=0; j<nbins; j++)
    {
        double x = -1.0 + (j+0.5)*2.0/nbins;
        chebyshev(x, n, l);
        double y = 0.0;
        for (unsigned int i = 0; i<n+1; i++){
            y += coeffs.at(i) * l.at(i);
        }
        interpol->SetBinContent(j+1, y);
    }
    interpol->SetLineColor(2);
    interpol->Draw("same l");

    canvas->Update();
    canvas->Print("chebyshev.eps", "eps");
    delete data;
    delete interpol;
    //delete pol1d;
    delete canvas;
}

void test_legendre_2d()
{
    gStyle->SetPalette(1);
    TRandom3 rnd;
    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200) ;
    canvas->cd();
    unsigned int nbins = 200;
    TH2D * data = new TH2D("data", "data", nbins, -1.0, 1.0, nbins, -1.0, 1.0);
    for (unsigned int i =0; i<1000000; i++)
    {
        //double r = rnd.Gaus(0.7,0.2);
        //double phi = rnd.Rndm() * 2.0 * TMath::Pi();
        //data->Fill(fabs(r)*cos(phi),fabs(r)*sin(phi));
        data->Fill(rnd.Gaus(-0.1,0.3), rnd.Gaus(-0.1,0.3));
        data->Fill(rnd.Gaus(0.5,0.3), rnd.Gaus(0.5,0.3));
    }
    //data->FillRandom("gaus", 200000);
    data->Draw("colz");
    canvas->Update();
    canvas->Print("2d.eps", "eps");
    unsigned int n = 10;
    std::vector<double> legendre_x, legendre_y;
    double dv = 2.0/nbins*2.0/nbins;
    std::vector< std::vector<double > > coeffs(n+1, std::vector<double >(n+1, 0.0));
    for (unsigned int i=0; i<nbins; i++)
    {
        double x = -1.0 + (i+0.5)*2.0/nbins;
        for (unsigned int j=0; j<nbins; j++)
        {
            double y = -1.0 + (j+0.5)*2.0/nbins;
            double z = data->GetBinContent(i+1,j+1);
            //orthonormal_legendre(x, n, legendre_x);
            //orthonormal_legendre(y, n, legendre_y);
            legendre(x, n, legendre_x);
            legendre(y, n, legendre_y);
            for (unsigned int l = 0; l<n+1; l++)
                for (unsigned int m = 0; m<n+1; m++)
                    coeffs.at(l).at(m) += z*legendre_x.at(l)*legendre_y.at(m)*dv* (0.5*(2.0*l+1.0)) * (0.5*(2.0*m+1.0));
        }
    }
    for (unsigned int l = 0; l<n+1; l++)
        for (unsigned int m = 0; m<n+1; m++)
            spdlog::debug("coefficient ({0:d}, {1:d}) = {2:d}", l, m, coeffs.at(l).at(m));


    TH2D * interpol = new TH2D("inter", "inter", nbins, -1.0, 1.0, nbins, -1.0, 1.0);
    for (unsigned int i=0; i<nbins; i++)
    {
        double x = -1.0 + (i+0.5)*2.0/nbins;
        for (unsigned int j=0; j<nbins; j++)
        {
            double y = -1.0 + (j+0.5)*2.0/nbins;
            //orthonormal_legendre(x, n, legendre_x);
            //orthonormal_legendre(y, n, legendre_y);
            legendre(x, n, legendre_x);
            legendre(y, n, legendre_y);

            double z = 0.0;
            for (unsigned int l = 0; l<n+1; l++)
                for (unsigned int m = 0; m<n+1; m++)
                    z += coeffs.at(l).at(m)*legendre_x.at(l)*legendre_y.at(m);
            interpol->SetBinContent(i+1, j+1, z);
        }
    }
    //calculate chi^2
    double chi2 = 0.0;
    for (unsigned int i=0; i<nbins; i++)
    {
        for (unsigned int j=0; j<nbins; j++)
        {
            double d = data->GetBinContent(i+1, j+1);
            double dd = data->GetBinError(i+1, j+1);
            double f = interpol->GetBinContent(i+1, j+1);
            if (dd != 0.0)
                chi2 += (d-f)*(d-f)/dd/dd;
        }
    }
    spdlog::debug("Found a chi^2 of {0:f with {1:d} bins.", chi2, nbins*nbins);
    //cleanup
    interpol->Draw("colz");
    canvas->Update();
    canvas->Print("2dinter.eps", "eps");
    delete data;
    delete interpol;
    delete canvas;
}

/** This function tests the different versions of the complex error function.
   * Unfortunately the complex error function is not yet in the standard c library.
   **/
void test_errf()
{
    gStyle->SetPalette(1);
    unsigned int divisions=1000;
    double min=-10.0;
    double max=+10.0;
    TH2D* real_1 = new TH2D("real_1", "real part", divisions, min, max, divisions, min, max);
    real_1->SetMinimum(-20.0);
    real_1->SetMaximum(+20.0);
    TH2D* imaginary_1 = new TH2D("imaginary_1", "imaginary part", divisions, min, max, divisions, min, max);
    imaginary_1->SetMinimum(-20.0);
    imaginary_1->SetMaximum(+20.0);
    TH2D* real_2 = new TH2D("real_2", "real part", divisions, min, max, divisions, min, max);
    TH2D* imaginary_2 = new TH2D("imaginary_2", "imaginary part", divisions, min, max, divisions, min, max);
    TH2D* real_3 = new TH2D("real_3", "real part", divisions, min, max, divisions, min, max);
    TH2D* imaginary_3 = new TH2D("imaginary_3", "imaginary part", divisions, min, max, divisions, min, max);
    double r, i;
    std::complex<double> result, arg;
    for (unsigned int re = 0; re < divisions; re++)
        for (unsigned int im = 0; im < divisions; im++)
        {
            r = min + (max - min)/divisions*re;
            i = min + (max - min)/divisions*im;
            arg = std::complex<double>(r, i);
            result = 1.0-cErrF(arg);
            //if (!std::isnan(result.real()))
            real_1->SetBinContent(re, im, result.real());
            imaginary_1->SetBinContent(re, im, result.imag());
            result = 1.0-cErrF_2(arg);
            //if (!std::isnan(result.real()))
            real_2->SetBinContent(re, im, result.real());
            imaginary_2->SetBinContent(re, im, result.imag());
            result = 1.0-cErrF_3(arg);
            //if (!std::isnan(result.real()))
            real_3->SetBinContent(re, im, result.real());
            imaginary_3->SetBinContent(re, im, result.imag());
        }

    TCanvas *canvas_1 = new TCanvas("canvas_1", "canvas 1",1600,1200) ;
    canvas_1->Divide(2,1);
    canvas_1->cd(1);
    real_1->Draw("COLZ");
    canvas_1->cd(2);
    imaginary_1->Draw("COLZ");
    canvas_1->Update();
    canvas_1->Print("error_1.eps", "eps");

    TCanvas *canvas_2 = new TCanvas("canvas_2", "canvas 2",1600,1200) ;
    canvas_2->Divide(2,1);
    canvas_2->cd(1);
    real_2->Draw("COLZ");
    canvas_2->cd(2);
    imaginary_2->Draw("COLZ");
    canvas_2->Update();
    canvas_2->Print("error_2.eps", "eps");

    TCanvas *canvas_3 = new TCanvas("canvas_3", "canvas 3",1600,1200) ;
    canvas_3->Divide(2,1);
    canvas_3->cd(1);
    real_3->Draw("COLZ");
    canvas_3->cd(2);
    imaginary_3->Draw("COLZ");
    canvas_3->Update();
    canvas_3->Print("error_3.eps", "eps");

    TH2D* real_12 = new TH2D((*real_1) / (*real_2));
    TH2D* imaginary_12 = new TH2D((*imaginary_1) / (*imaginary_2));
    TCanvas *canvas_12 = new TCanvas("canvas_12", "canvas 1 - 2",1600,1200) ;
    canvas_12->Divide(2,1);
    canvas_12->cd(1);
    real_12->Draw("COLZ");
    canvas_12->cd(2);
    imaginary_12->Draw("COLZ");
    canvas_12->Update();
    canvas_12->Print("error_1_over_2.eps", "eps");

    TH2D* real_23 = new TH2D((*real_2) / (*real_3));
    TH2D* imaginary_23 = new TH2D((*imaginary_2) / (*imaginary_3));
    TCanvas *canvas_23 = new TCanvas("canvas_23", "canvas 1 - 2",1600,1200) ;
    canvas_23->Divide(2,1);
    canvas_23->cd(1);
    real_23->Draw("COLZ");
    canvas_23->cd(2);
    imaginary_23->Draw("COLZ");
    canvas_23->Update();
    canvas_23->Print("error_2_over_3.eps", "eps");

    TH2D* real_13 = new TH2D((*real_1) / (*real_3));
    TH2D* imaginary_13 = new TH2D((*imaginary_1) / (*imaginary_3));
    TCanvas *canvas_13 = new TCanvas("canvas_13", "canvas 1 - 2",1600,1200) ;
    canvas_13->Divide(2,1);
    canvas_13->cd(1);
    real_13->Draw("COLZ");
    canvas_13->cd(2);
    imaginary_13->Draw("COLZ");
    canvas_13->Update();
    canvas_13->Print("error_1_over_3.eps", "eps");

}


void test_convolutions()
{
    gStyle->SetPalette(1);
    gStyle->SetOptFit(0);
    gStyle->SetOptStat(0);

    unsigned int divisions=1000;
    double min=-TMath::Pi();
    double max=+2.0*TMath::Pi();
    TH1D* conv_sin = new TH1D("conv_sin", "", divisions, min, max);
    TH1D* conv_cos = new TH1D("conv_cos", "", divisions, min, max);
    TH1D* conv_sin_2 = new TH1D("conv_sin_2", "", divisions, min, max);
    TH1D* conv_cos_2 = new TH1D("conv_cos_2", "", divisions, min, max);
    TH1D* conv_2_sin = new TH1D("conv_2_sin", "", divisions, min, max);
    TH1D* conv_2_cos = new TH1D("conv_2_cos", "", divisions, min, max);

    TH1D* conv_cos_2_sin = new TH1D("conv_cos_2_sin", "", divisions, min, max);
    TH1D* conv_sin_3 = new TH1D("conv_sin_3", "", divisions, min, max);
    TH1D* conv_2_sin_sin = new TH1D("conv_2_sin_sin", "", divisions, min, max);

    double sigma = 0.03;
    for (unsigned int i = 0; i < divisions; i++)
    {
        double angle = min+(max - min)/double(divisions)*i;
        conv_sin->SetBinContent(i+1, convoluted_sin(angle, sigma));
        conv_cos->SetBinContent(i+1, convoluted_cos(angle, sigma));
        conv_sin_2->SetBinContent(i+1, convoluted_sin_2(angle, sigma));
        conv_cos_2->SetBinContent(i+1, convoluted_cos_2(angle, sigma));
        conv_2_sin->SetBinContent(i+1, convoluted_2_sin(angle, sigma));
        conv_2_cos->SetBinContent(i+1, convoluted_2_cos(angle, sigma));

        conv_cos_2_sin->SetBinContent(i+1, convoluted_cos_2_sin(angle, sigma));
        conv_sin_3->SetBinContent(i+1, convoluted_sin_3(angle, sigma));
        //std::cout << convoluted_sin_3(angle, sigma) << std::endl;
        conv_2_sin_sin->SetBinContent(i+1, convoluted_2_sin_sin(angle, sigma));

    }

    conv_sin->SetLineColor(2);
    conv_cos->SetLineColor(3);
    conv_sin_2->SetLineColor(4);
    conv_cos_2->SetLineColor(5);
    conv_2_sin->SetLineColor(6);
    conv_2_cos->SetLineColor(7);
    conv_cos_2_sin->SetLineColor(8);
    conv_sin_3->SetLineColor(9);
    conv_2_sin_sin->SetLineColor(15);

    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200) ;
    canvas->cd();
    conv_sin->SetMinimum(-1.0);
    conv_sin->SetMaximum(1.0);
    conv_sin->Draw();
    conv_cos->Draw("same");
    conv_sin_2->Draw("same");
    conv_cos_2->Draw("same");
    conv_2_sin->Draw("same");
    conv_2_cos->Draw("same");
    conv_cos_2_sin->Draw("same");
    conv_sin_3->Draw("same");
    conv_2_sin_sin->Draw("same");
    TLegend* conv_legend = new TLegend(0.7,0.6,0.9,0.9);
    conv_legend->AddEntry(conv_sin, "sin #Theta", "L");
    conv_legend->AddEntry(conv_cos, "cos #Theta", "L");
    conv_legend->AddEntry(conv_sin_2, "sin #Theta^{2}", "L");
    conv_legend->AddEntry(conv_cos_2, "cos #Theta^{2}", "L");
    conv_legend->AddEntry(conv_2_sin, "sin 2#Theta", "L");
    conv_legend->AddEntry(conv_2_cos, "cos 2#Theta", "L");
    conv_legend->AddEntry(conv_cos_2_sin, "cos #Theta^{2} sin #Theta", "L");
    conv_legend->AddEntry(conv_sin_3, "sin #Theta ^{3}", "L");
    conv_legend->AddEntry(conv_2_sin_sin, "sin 2#Theta sin #Theta", "L");
    conv_legend->Draw();
    canvas->Print("convolutions.eps", "eps");

    delete conv_sin;
    delete conv_cos;
    delete conv_sin_2;
    delete conv_cos_2;
    delete conv_2_sin;
    delete conv_2_cos;
    delete conv_cos_2_sin;
    delete conv_sin_3;
    delete conv_2_sin_sin;

    delete canvas;

    TH1D* conv_cos_2_a = new TH1D("conv_cos_2_a", "", divisions, -1.0, +1.0);
    TH1D* conv_cos_2_b = new TH1D("conv_cos_2_b", "", divisions, -1.0, +1.0);
    TH1D* conv_cos_2_c = new TH1D("conv_cos_2_c", "", divisions, -1.0, +1.0);
    TH1D* conv_cos_2_d = new TH1D("conv_cos_2_d", "", divisions, -1.0, +1.0);
    TH1D* conv_cos_2_e = new TH1D("conv_cos_2_e", "", divisions, -1.0, +1.0);
    conv_cos_2_a->SetMinimum(-1.0);
    conv_cos_2_a->SetMaximum(1.0);
    for (unsigned int i = 0; i < divisions; i++)
    {
        double cos_theta = -1.0 + 2.0/divisions*(i+0.5);
        double angle = acos(cos_theta);
        conv_cos_2_a->SetBinContent(i+1, cos_theta*cos_theta);
        conv_cos_2_b->SetBinContent(i+1, convoluted_cos_2(angle, sigma));
        conv_cos_2_c->SetBinContent(i+1, convoluted_cos_2(angle, sigma)
                                    +convoluted_cos_2(-angle, sigma)
                                    +convoluted_cos_2(2.0*TMath::Pi()-angle, sigma));
        conv_cos_2_d->SetBinContent(i+1,1.0/sin(angle)*convoluted_cos_2_sin(angle, sigma));
        conv_cos_2_e->SetBinContent(i+1,1.0/sin(angle)*(convoluted_cos_2_sin(angle, sigma)
                                                        +convoluted_cos_2_sin(-angle, sigma)
                                                        +convoluted_cos_2_sin(2.0*TMath::Pi()-angle, sigma)
                                                        ));
    }

    TCanvas *canvas2 = new TCanvas("canvas2", "canvas",1600,1200);
    canvas2->cd();
    conv_cos_2_b->SetLineColor(2);
    conv_cos_2_c->SetLineColor(3);
    conv_cos_2_d->SetLineColor(4);
    conv_cos_2_e->SetLineColor(5);
    conv_cos_2_a->Draw();
    conv_cos_2_b->Draw("same");
    conv_cos_2_c->Draw("same");
    conv_cos_2_d->Draw("same");
    conv_cos_2_e->Draw("same");
    canvas2->Print("overcos.eps", "eps");
    delete conv_cos_2_a;
    delete conv_cos_2_b;
    delete conv_cos_2_c;
    delete conv_cos_2_d;
    delete conv_cos_2_e;
    delete canvas2;
}

void test_int_convolutions()
{
    gStyle->SetPalette(1);
    gStyle->SetOptFit(0);
    gStyle->SetOptStat(0);

    unsigned int divisions=1000;
    double min=-TMath::Pi();
    double max=+2.0*TMath::Pi();
    TH1D* conv_sin = new TH1D("conv_sin", "", divisions, min, max);
    TH1D* conv_cos = new TH1D("conv_cos", "", divisions, min, max);
    TH1D* conv_sin_2 = new TH1D("conv_sin_2", "", divisions, min, max);
    TH1D* conv_cos_2 = new TH1D("conv_cos_2", "", divisions, min, max);
    TH1D* conv_2_sin = new TH1D("conv_2_sin", "", divisions, min, max);
    TH1D* conv_2_cos = new TH1D("conv_2_cos", "", divisions, min, max);

    TH1D* conv_cos_2_sin = new TH1D("conv_cos_2_sin", "", divisions, min, max);
    TH1D* conv_sin_3 = new TH1D("conv_sin_3", "", divisions, min, max);
    TH1D* conv_2_sin_sin = new TH1D("conv_2_sin_sin", "", divisions, min, max);

    double sigma = 0.03;
    for (unsigned int i = 0; i < divisions; i++)
    {
        double angle = min+(max - min)/double(divisions)*i;
        conv_sin->SetBinContent(i+1, int_convoluted_sin(angle, sigma)-int_convoluted_sin(min, sigma));

        // conv_cos->SetBinContent(i+1, convoluted_cos(angle, sigma));
        conv_sin_2->SetBinContent(i+1, int_convoluted_sin_2(angle, sigma) - int_convoluted_sin_2(min, sigma));
        // conv_cos_2->SetBinContent(i+1, convoluted_cos_2(angle, sigma));
        conv_2_sin->SetBinContent(i+1, int_convoluted_2_sin(angle, sigma)-int_convoluted_2_sin(min, sigma));
        // conv_2_cos->SetBinContent(i+1, convoluted_2_cos(angle, sigma));

        conv_cos_2_sin->SetBinContent(i+1, int_convoluted_cos_2_sin(angle, sigma) - int_convoluted_cos_2_sin(min, sigma));
        conv_sin_3->SetBinContent(i+1, int_convoluted_sin_3(angle, sigma)-int_convoluted_sin_3(min, sigma));
        conv_2_sin_sin->SetBinContent(i+1, int_convoluted_2_sin_sin(angle, sigma) - int_convoluted_2_sin_sin(min, sigma));
    }

    conv_sin->SetLineColor(2);
    conv_cos->SetLineColor(3);
    conv_sin_2->SetLineColor(4);
    conv_cos_2->SetLineColor(5);
    conv_2_sin->SetLineColor(6);
    conv_2_cos->SetLineColor(7);
    conv_cos_2_sin->SetLineColor(8);
    conv_sin_3->SetLineColor(9);
    conv_2_sin_sin->SetLineColor(15);

    TCanvas *canvas = new TCanvas("canvas", "canvas",1600,1200) ;
    canvas->cd();
    conv_sin->SetMinimum(-0.5);
    conv_sin->SetMaximum(3.0);
    conv_sin->Draw();

    // conv_cos->Draw("same");
    conv_sin_2->Draw("same");
    // conv_cos_2->Draw("same");
    conv_2_sin->Draw("same");
    // conv_2_cos->Draw("same");
    conv_cos_2_sin->Draw("same");
    conv_sin_3->Draw("same");
    conv_2_sin_sin->Draw("same");
    TLegend* conv_legend = new TLegend(0.7,0.6,0.9,0.9);
    conv_legend->AddEntry(conv_sin, "sin #Theta", "L");
    conv_legend->AddEntry(conv_cos, "cos #Theta", "L");
    conv_legend->AddEntry(conv_sin_2, "sin #Theta^{2}", "L");
    conv_legend->AddEntry(conv_cos_2, "cos #Theta^{2}", "L");
    conv_legend->AddEntry(conv_2_sin, "sin 2#Theta", "L");
    conv_legend->AddEntry(conv_2_cos, "cos 2#Theta", "L");
    conv_legend->AddEntry(conv_cos_2_sin, "cos #Theta^{2} sin #Theta", "L");
    conv_legend->AddEntry(conv_sin_3, "sin #Theta ^{3}", "L");
    conv_legend->AddEntry(conv_2_sin_sin, "sin 2#Theta sin #Theta", "L");
    conv_legend->Draw();
    canvas->Print("int_convolutions.eps", "eps");

    delete conv_sin;
    delete conv_cos;
    delete conv_sin_2;
    delete conv_cos_2;
    delete conv_2_sin;
    delete conv_2_cos;
    delete conv_cos_2_sin;
    delete conv_sin_3;
    delete conv_2_sin_sin;
    delete canvas;
}


}