Feature/fitting (#1)

Add root fitting code
Add adjusted and normalised datasets
Add mathematica fitting code
Add error estimate
Add chi square calculator

Pipfile

@@ -1,2 +1,6 @@
 [scripts]
 "build:report" = "sh ./scripts/build:report.sh"
+
+[packages]
+jupyter = "*"
+metakernel = "*"

code/README.md

@@ -0,0 +1,3 @@
+# Code
+This folder contains all the code we used.
+Still a work in progress.

code/fitting/ROOT/functions.hpp

@@ -0,0 +1,59 @@
+#ifndef FUNCTIONS_HPP
+#define FUNCTIONS_HPP
+
+inline TCanvas* cv;
+inline TF1* r_pi_fun;
+inline TF1* r_sigma_fun;
+inline TF1* R_pi_fun;
+inline TF1* R_sigma_fun;
+inline TF1* R_pi_raw_fun;
+inline TF1* R_sigma_raw_fun;
+
+// Ampiezza riflessa con polarizzazione pi
+inline double r_pi(double *xx,double *par) {
+  const double n1         =  par[0];
+  const double n2         =  par[1];
+  const double deltaTheta =  par[2];
+
+  const double theta1 = (xx[0] - deltaTheta) / 180 * TMath::Pi();
+
+  const double theta2 = TMath::ASin(n1 / n2 * TMath::Sin(theta1));
+
+  return (n1 * cos(theta2) - n2 * cos(theta1)) / (n1 * cos(theta2) + n2 * cos(theta1));
+}
+
+// Ampiezza riflessa con polarizzazione sigma
+inline double r_sigma(double *xx,double *par) {
+  const double n1         =  par[0];
+  const double n2         =  par[1];
+  const double deltaTheta =  par[2];
+
+  const double theta1 = (xx[0] - deltaTheta) / 180 * TMath::Pi();
+
+  const double theta2 = TMath::ASin(n1 / n2 * TMath::Sin(theta1));
+
+  return (n1 * cos(theta1) - n2 * cos(theta2)) / (n1 * cos(theta1) + n2 * cos(theta2));
+}
+
+// Intensità riflessa con polarizzazione pi
+inline double R_pi(double *xx,double *par) {
+  const double I_0 = par[3];
+
+  return I_0 * pow(r_pi(xx, par), 2);
+}
+
+// Ampiezza riflessa con polarizzazione sigma
+inline double R_sigma(double *xx,double *par) {
+  const double I_0 = par[3];
+
+  return I_0 * pow(r_sigma(xx, par), 2);
+}
+
+//double T_p2(double *xx,double *par) {
+//  return 1 - R_p2(xx, par);
+//}
+//
+//double T_s2(double *xx,double *par) {
+//  return 1 - R_s2(xx, par);
+
+#endif // define FUNCTIONS_HPP

code/fitting/ROOT/graphs.hpp

@@ -0,0 +1,7 @@
+#ifndef GRAPHS_HPP
+#define GRAPHS_HPP
+
+inline TGraphErrors* R_pi_raw_graph;
+inline TGraphErrors* R_sigma_raw_graph;
+
+#endif // define GRAPHS_HPP

code/fitting/ROOT/macro.cpp

@@ -0,0 +1,62 @@
+#include "functions.hpp"
+#include "graphs.hpp"
+
+void macro() {
+  gStyle->SetOptFit(111111);
+
+  cv = new TCanvas("rawdata", "Raw Data");
+  cv->Divide(2, 2);
+
+  cv->cd(1);
+  R_pi_raw_graph = new TGraphErrors("../../data/pi/run3-adjusted.csv", "%lg %lg %lg %lg", "  ");
+  R_pi_raw_graph->Draw();
+
+  R_pi_raw_fun = new TF1("R_pi_raw_fun", R_pi, 0, 90, 4);
+  R_pi_raw_fun->SetParNames("$n_1$", "$n_2$", "$$\\Delta\\theta$$", "$$I_0$$");
+  R_pi_raw_fun->SetParameters(1, 1.519, 0, 1000);
+  R_pi_raw_fun->SetParLimits(0, 1, 1);
+  R_pi_raw_fun->SetParLimits(2, -0.00001, 0.00001);
+  R_pi_raw_fun->Draw("same");
+  R_pi_raw_graph->Fit(R_pi_raw_fun);
+
+  cv->cd(3);
+  ifstream ifs("../../data/pi/run3-adjusted.csv");
+  int npoints = 38;
+  double x[npoints];
+  double y[npoints];
+  for (int i = 0; i < npoints; ++i) {
+    double _;
+    ifs >> x[i] >> y[i] >> _ >> _;
+    y[i] -= R_pi_raw_fun->Eval(x[i]);
+  }
+  auto pi_residual = new TGraph(npoints, x, y);
+  pi_residual->Draw();
+
+  // Sigma polarisation ////////////////////////////////////////////////////////
+  cv->cd(2);
+  R_sigma_raw_graph = new TGraphErrors("../../data/sigma/run1-adjusted.csv", "%lg %lg %lg %lg");
+  R_sigma_raw_graph->Draw();
+
+  R_sigma_raw_fun = new TF1("R_sigma_raw_fun", R_sigma, 0, 90, 4);
+  R_sigma_raw_fun->SetParNames("$n_1$", "$n_2$", "$$\\Delta\\theta$$", "$$I_0$$");
+  R_sigma_raw_fun->SetParameters(1, 1.519, 0, 1000);
+  R_sigma_raw_fun->SetParLimits(0, 1, 1);
+  R_sigma_raw_fun->SetParLimits(2, -0.00001, 0.00001);
+  R_sigma_raw_fun->Draw("same");
+  R_sigma_raw_graph->Fit(R_sigma_raw_fun);
+
+  cv->cd(4);
+  ifs.close();
+  ifs.open("../../data/sigma/run1-adjusted.csv");
+  for (int i = 0; i < npoints; ++i) {
+    double _;
+    ifs >> x[i] >> y[i] >> _ >> _;
+    y[i] -= R_sigma_raw_fun->Eval(x[i]);
+  }
+  auto sigma_resiudal = new TGraph(npoints, x, y);
+  sigma_resiudal->Draw();
+}
+
+
+
+

code/fitting/chi-square.py

@@ -0,0 +1,55 @@
+from math import pi, cos, asin, sin
+
+n2_pi = 1.493
+n2_sigma = 1.490
+
+I_0pi = 2550
+I_0sigma = 1340
+
+a = 0.29
+theta_brewster = 56.8
+
+def R_pi(x):
+  theta1 = x / 180 * pi
+  theta2 = asin(1 / n2_pi * sin(theta1))
+
+  r_pi = (1 * cos(theta2) - n2_pi * cos(theta1)) / (1 * cos(theta2) + n2_pi * cos(theta1))
+  return I_0pi * (r_pi ** 2)
+
+def R_sigma(x):
+  theta1 = x / 180 * pi
+  theta2 = asin(1 / n2_pi * sin(theta1))
+
+  r_sigma = (1 * cos(theta1) - n2_sigma * cos(theta2)) / (1 * cos(theta1) + n2_sigma * cos(theta2))
+  return I_0sigma * (r_sigma ** 2)
+
+def D(func, x):
+  return (func(x+0.0000001) - func(x-0.0000001)) / 0.0000002
+
+def parabola(x):
+  return a*((x - theta_brewster)**2)
+
+def calc_chisquare(xs, ys, delta_xs, delta_ys, func):
+  delta_y1s = [delta_xs * D(func, x) for x, delta_xs in zip(xs, delta_xs)]
+  delta_ys_corrected = [dy + dy1 for dy, dy1 in zip(delta_ys, delta_y1s)]
+  os = [func(x) for x in xs]
+  chi_square = sum([((e - o) / dyc)**2 for e, o, dyc in zip(ys, os, delta_ys_corrected)])
+  return chi_square
+
+
+datasets = ["./data/pi/run3-adjusted.csv", "./data/sigma/run1-adjusted.csv", "./data/pi/run3-parabola.csv"]
+funcs = [R_pi, R_sigma, parabola]
+constraints = [2, 2, 2]
+for dataset, func, constraints in zip(datasets, funcs, constraints):
+  xs = []
+  ys = []
+  delta_xs = []
+  delta_ys = []
+  with open(dataset) as file:
+    for line in file:
+      x, y, dx, dy = line.split()
+      xs.append(float(x))
+      ys.append(float(y))
+      delta_xs.append(0) # float(dx)) we don't want to recalculate this.
+      delta_ys.append(float(dy))
+  print(dataset, calc_chisquare(xs, ys, delta_xs, delta_ys, func) / (len(xs) - constraints))