kopecna
/
EWP-BplusToKstMuMu-AngAna
2
1
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Angular analysis of B+->K*+(K+pi0)mumu
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
29 Commits
2 Branches
0 Tags
92 MiB
 
 
 
 
Branch: master
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from 'master'
${ noResults }
EWP-BplusToKstMuMu-AngAna/Code/FCNCFitter/sources/tests.cc

781 lines
26 KiB
Raw Permalink Blame History

//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;
}
}