Esempi: Appendice 4
Contents
Esempi: Appendice 4¶
A4.0 generazione di coppie di punti secondo una legge funzionale y = g(x), con dispersione gaussiana lungo y¶
main_00.cpp
/*
c++ -o main_00 main_00.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "generazione.h"
using namespace std ;
// retta con termine noto non nullo
double g (double x)
{
return 3.14 + 2 * x ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv)
{
if (argc < 3)
{
cout << "usage: " << argv[0] << " sigma_y nomefile.txt [Npoints = 5]" << endl ;
exit (1) ;
}
double sigma = fabs (atof (argv[1])) ;
ofstream output_file ;
output_file.open (argv[2]) ;
int Npoints = 5 ;
if (argc > 3) Npoints = atoi (argv[3]) ;
for (int i = 0 ; i < Npoints ; ++i)
{
double epsilon = rand_TAC_gaus (sigma) ;
output_file << i << " " << g (i) + epsilon << "\n" ;
}
output_file.close () ;
return 0 ;
}
A4.1 determinazione dei parametri di una funzione lineare y = a + bx con il metodo dei minimi quadrati¶
main_01.cpp
/*
c++ -o main_01 algebra_2.cc main_01.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
#include "algebra_2.h"
using namespace std ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv)
{
if (argc < 3)
{
cout << "usage: " << argv[0] << " sigma_y nomefile.txt" << endl ;
exit (1) ;
}
// lettura del file di input
// -------------------------
double sigma = atof (argv[1]) ;
ifstream input_file ;
input_file.open (argv[2]) ;
vector<double> asse_x ;
vector<double> asse_y ;
while (true)
{
double dummy_x = 0. ;
double dummy_y = 0. ;
input_file >> dummy_x ;
input_file >> dummy_y ;
if (input_file.eof () == true) break ;
asse_x.push_back (dummy_x) ;
asse_y.push_back (dummy_y) ;
}
input_file.close () ;
int N_points = asse_x.size () ;
// creazione delle matrici del metodo dei minimi quadrati
// -------------------------
matrice H (N_points, 2) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
H.setCoord (i_point, 0, 1) ;
H.setCoord (i_point, 1, asse_x.at (i_point)) ;
}
vettore y (asse_y) ;
// ipotesi: incertezza costante su tutti i valori
matrice V (N_points) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
V.setCoord (i_point, i_point, sigma * sigma) ;
// calcolo dei parametri di interesse
// -------------------------
matrice V_inv = V.inversa () ;
//Metzger, cap. 8.5.2, eq. 8.116
matrice theta_v = (H.trasposta () * V_inv * H).inversa () ;
//Metzger, cap. 8.5.2, eq. 8.112
vettore theta = (theta_v * (H.trasposta () * V_inv)) * y ;
cout << "parametri risultanti dai mimimi quadrati: \n" ;
theta.stampa () ;
cout << "matrice di covarianza dei parametri risultanti dai mimimi quadrati: \n" ;
theta_v.stampa () ;
cout << "termine noto: " << theta.at (0) << " +- " << sqrt (theta_v.at (0, 0)) << endl ;
cout << "pendenza: " << theta.at (1) << " +- " << sqrt (theta_v.at (1, 1)) << endl ;
return 0 ;
}
algebra_2.h
#ifndef algebra_2_h
#define algebra_2_h
#include <cmath>
#include <vector>
class matrice ;
class vettore
{
public:
vettore (int N) ;
vettore (const std::vector<double> & v) ;
vettore (const vettore & orig) ;
vettore & operator = (const vettore & orig) ;
~vettore () ;
void setCoord (int i, double val) ;
double norm () const ;
int N () const ;
double at (int i) const ;
void stampa () const ;
double operator[] (int i) const ;
vettore operator+ (const vettore & v) const ;
vettore operator- (const vettore & v) const ;
vettore operator* (double val) const ;
double dot (const vettore & v) const ;
private:
double * m_elementi ;
int m_N ;
} ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
class matrice
{
public:
matrice (int R) ;
matrice (int R, int C) ;
matrice (const matrice & orig) ;
matrice & operator= (const matrice & orig) ;
~matrice () ;
void setCoord (int i, int j, double val) ;
double at (int i, int j) const ;
void stampa () const ;
bool quadrata () const ;
int rango () const ;
int N_righe () const ;
int N_colonne () const ;
bool simmetrica () const ;
void dimensioni () const ;
matrice minore (int r, int c) const ;
matrice inversa () const ;
matrice trasposta () const ;
double determinante () const ;
void operator*= (double val) ;
private:
int index (int i, int j) const ;
int m_R ;
int m_C ;
double * m_elementi ;
} ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
vettore operator* (const matrice & M, const vettore & v) ;
matrice operator* (const matrice & M1, const matrice & M2) ;
matrice operator+ (const matrice & M1, const matrice & M2) ;
#endif
algebra_2.cc
#include "algebra_2.h"
#include <iostream>
using namespace std ;
vettore::vettore (int N) :
m_N (N),
m_elementi (new double[N])
{
for (int i = 0 ; i < m_N ; ++i) m_elementi[i] = 0. ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore::vettore (const std::vector<double> & v) :
m_N (v.size ()),
m_elementi (new double[v.size ()])
{
for (int i = 0 ; i < m_N ; ++i) m_elementi[i] = v.at (i) ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore::vettore (const vettore & orig) :
m_N (orig.m_N),
m_elementi (new double[orig.m_N])
{
for (int i = 0 ; i < m_N ; ++i) m_elementi[i] = orig.m_elementi[i] ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore & vettore::operator = (const vettore & orig)
{
m_N = orig.m_N ;
m_elementi = new double[m_N] ;
for (int i = 0 ; i < m_N ; ++i) m_elementi[i] = orig.m_elementi[i] ;
return *this ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore::~vettore () { delete [] m_elementi ; }
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void vettore::setCoord (int i, double val)
{
if (i < 0 || i > m_N-1) return ;
m_elementi[i] = val ;
return ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double vettore::norm () const
{
double sum = m_elementi[0] * m_elementi[0] ;
for (int i = 1 ; i < m_N ; ++i) sum += m_elementi[i] * m_elementi[i] ;
return sqrt (sum) ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
int vettore::N () const {return m_N ;}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double vettore::at (int i) const
{
if (i < 0 || i > m_N-1) return 0. ;
return m_elementi[i] ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void vettore::stampa () const
{
for (int i = 0 ; i < m_N ; ++i)
cout << this->at (i) << "\t" ;
cout << "\n" ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double vettore::operator[] (int i) const
{
return m_elementi[i] ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore vettore::operator+ (const vettore & v) const
{
vettore result (*this) ;
for (int i = 0 ; i < m_N ; ++i)
result.setCoord (i, result.at (i) + v.at (i)) ;
return result ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore vettore::operator- (const vettore & v) const
{
vettore result (*this) ;
for (int i = 0 ; i < m_N ; ++i)
result.setCoord (i, result.at (i) - v.at (i)) ;
return result ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
vettore vettore::operator* (double val) const
{
vettore result (*this) ;
for (int i = 0 ; i < m_N ; ++i)
result.setCoord (i, result.at (i) * val) ;
return result ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double vettore::dot (const vettore & v) const
{
if (m_N != v.N ())
{
cerr << "prodotto scalare fra vettori di dimensione diversa" << endl ;
exit (1) ;
}
double result = 0 ;
for (int i = 0 ; i < m_N ; ++i) result += this->at (i) * v.at (i) ;
return result ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
matrice::matrice (int R) :
m_R (R),
m_C (R),
m_elementi (new double [R * R])
{
for (int i = 0 ; i < m_R ; ++i)
for (int j = 0 ; j < m_C ; ++j)
this->setCoord (i, j, 0.) ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice::matrice (int R, int C) :
m_R (R),
m_C (C),
m_elementi (new double [R * C])
{
for (int i = 0 ; i < m_R ; ++i)
for (int j = 0 ; j < m_C ; ++j)
this->setCoord (i, j, 0.) ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice::matrice (const matrice & orig) :
m_R (orig.m_R),
m_C (orig.m_C),
m_elementi (new double [orig.m_R * orig.m_C])
{
for (int i = 0 ; i < m_R ; ++i)
for (int j = 0 ; j < m_C ; ++j)
this->setCoord (i, j, orig.at (i,j)) ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice & matrice::operator= (const matrice & orig)
{
m_R = orig.m_R ;
m_C = orig.m_C ;
m_elementi = new double [orig.m_R * orig.m_C] ;
for (int i = 0 ; i < m_R ; ++i)
for (int j = 0 ; j < m_C ; ++j)
this->setCoord (i, j, orig.at (i,j)) ;
return *this ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice::~matrice () { delete [] m_elementi ; }
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
int matrice::index (int i, int j) const
{
int index = j + i * m_C ;
if (index < 0 || index >= m_R * m_C)
{
cerr << "indici matrice fuori intervallo" << endl ;
exit (1) ;
}
return index ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void matrice::setCoord (int i, int j, double val)
{
m_elementi[index (i, j)] = val ;
return ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double matrice::at (int i, int j) const
{
return m_elementi[index (i, j)] ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void matrice::stampa () const
{
for (int i = 0 ; i < m_R ; ++i)
{
for (int j = 0 ; j < m_C ; ++j)
cout << this->at (i, j) << "\t" ;
cout << "\n" ;
}
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
bool matrice::quadrata () const {return (m_R == m_C) ;}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
// per matrici quadrate
int matrice::rango () const {return m_R ;}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
// per matrici quadrate
int matrice::N_righe () const {return m_R ;}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
// per matrici quadrate
int matrice::N_colonne () const {return m_C ;}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
bool matrice::simmetrica () const
{
if (!this->quadrata ()) return false ;
for (int i = 0 ; i < m_R ; ++i)
for (int j = i+1 ; j < m_R ; ++j)
if (fabs (this->at (i,j) - this->at (j,i)) > 0.00001) return false ;
return true ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void matrice::dimensioni () const { cout << m_R << " x " << m_C << endl ; }
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
// produce la sottomatrice in cui mancano riga r e colonna c
matrice matrice::minore (int r, int c) const
{
matrice M_out (m_R-1, m_C-1) ;
int i_new = 0 ;
for (int i = 0 ; i < m_R ; ++i)
{
if (i == r) continue ;
int j_new = 0 ;
for (int j = 0 ; j < m_C ; ++j)
{
if (j == c) continue ;
M_out.setCoord (i_new, j_new, this->at (i,j)) ;
++j_new ;
}
++i_new ;
}
return M_out ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice matrice::inversa () const
{
double invdet = this->determinante () ;
matrice M_out (m_C, m_R) ;
if (invdet < 0.000001) return M_out ;
invdet = 1. / invdet ;
for (int i = 0 ; i < m_R ; ++i)
{
for (int j = 0 ; j < m_C ; ++j)
{
M_out.setCoord (
j, i, // trasposta
this->minore (i,j).determinante () * pow (-1, i+j+2) // complemento algebrico
* invdet // divisione per il determinante
) ;
}
}
return M_out ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
matrice matrice::trasposta () const
{
matrice M_out (m_C, m_R) ;
for (int i = 0 ; i < m_R ; ++i)
for (int j = 0 ; j < m_C ; ++j)
M_out.setCoord (j, i, this->at (i,j)) ;
return M_out ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
double matrice::determinante () const
{
if (!this->quadrata ()) return 0. ;
if (this->rango () == 1) return this->at (0,0) ;
double det = 0. ;
for (int i = 0 ; i < m_R ; ++i)
{
det += this->at (0,i) * pow (-1, i + 2) * this->minore (0,i).determinante () ;
}
return det ;
}
// .... .... .... .... .... .... .... .... .... .... .... .... .... ....
void matrice::operator*= (double val)
{
for (int i = 0 ; i < m_R * m_C ; ++i)
{
m_elementi[i] *= val ;
}
return ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
vettore operator* (const matrice & M, const vettore & v)
{
if (M.N_colonne () != v.N ())
{
cerr << "prodotto fra matrice e vettore con dimensioni non compatibili" << endl ;
exit (1) ;
}
vettore w (M.N_righe ()) ;
for (int i = 0 ; i < M.N_righe () ; ++i)
{
double res = 0 ;
for (int j = 0 ; j < M.N_colonne () ; ++j) res += M.at (i,j) * v.at (j) ;
w.setCoord (i, res) ;
}
return w ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
matrice operator* (const matrice & M1, const matrice & M2)
{
if (M1.N_colonne () != M2.N_righe ())
{
cerr << "prodotto fra matrici con dimensioni non compatibili" << endl ;
exit (1) ;
}
matrice M_out (M1.N_righe (), M2.N_colonne ()) ;
for (int i = 0 ; i < M1.N_righe () ; ++i)
for (int k = 0 ; k < M2.N_colonne () ; ++k)
{
double res = 0 ;
for (int j = 0 ; j < M1.N_colonne () ; ++j) res += M1.at (i,j) * M2.at (j,k) ;
M_out.setCoord (i, k, res) ;
}
return M_out ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
matrice operator+ (const matrice & M1, const matrice & M2)
{
if (M1.N_colonne () != M2.N_colonne () ||
M1.N_righe () != M2.N_righe ())
{
cerr << "somma fra matrici con dimensioni non compatibili" << endl ;
exit (1) ;
}
matrice M_out (M1.N_righe (), M1.N_colonne ()) ;
for (int i = 0 ; i < M1.N_righe () ; ++i)
for (int k = 0 ; k < M1.N_colonne () ; ++k)
{
M_out.setCoord (i, k, M1.at (i,k) + M2.at (i,k)) ;
}
return M_out ;
}
generazione.h
#ifndef generazione_h
#define generazione_h
// media sempre zero
// non normalizzata, cosi' il massimo in zero vale sempre 1
double fgaus (double x, double sigma)
{
return exp (-0.5 * x * x / (sigma * sigma)) ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double rand_range (double min, double max)
{
return min + (max - min) * rand () / static_cast<double> (RAND_MAX) ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double rand_TAC_gaus (double sigma)
{
double x = 0. ;
double y = 0. ;
do {
x = rand_range (-5 * sigma, 5 * sigma) ;
y = rand_range (0, 1.) ;
} while (y > fgaus (x, sigma)) ;
return x ;
}
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
double rand_TAC_uniform (double sigma)
{
// sqrt (3) = 1.7320508076
return rand_range (-1.7320508076 * sigma, 1.7320508076 * sigma) ;
}
#endif
A4.2 studio del risultato del fit con toy experiment¶
main_02.cpp
/*
c++ -o main_02 `root-config --glibs --cflags` algebra_2.cc main_02.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
#include "TH1F.h"
#include "TH2F.h"
#include "TCanvas.h"
#include "algebra_2.h"
#include "generazione.h"
// retta con termine noto non nullo
double g (double x)
{
return 3.14 + 2 * x ;
}
using namespace std ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv)
{
if (argc < 4)
{
cout << "usage: " << argv[0] << " sigma_y N_points Ntoys" << endl ;
exit (1) ;
}
double sigma = atof (argv[1]) ;
int N_points = atoi (argv[2]) ;
int N_toys = atoi (argv[3]) ;
// istogrammi per il disegno dei risultati del fit
TH1F h_a ("h_a", "termine noto",
100, 3.14 * (1. - 1. * sigma), 3.14 * (1. + 1. * sigma) ) ;
TH1F h_b ("h_b", "coefficiente angolare",
100, 2. * (1. - 1. * sigma), 2. * (1. + 1. * sigma) ) ;
TH2F h_ab ("h_ab", "parametri",
50, 3.14 * (1. - 1. * sigma), 3.14 * (1. + 1. * sigma),
50, 2. * (1. - 1. * sigma), 2. * (1. + 1. * sigma) ) ;
h_ab.GetXaxis ()->SetTitle ("termine noto") ;
h_ab.GetYaxis ()->SetTitle ("coefficiente angolare") ;
h_ab.SetStats (0) ;
// contatori per la determinazione dell'intervallo di copertura
int cont_a = 0 ;
int cont_b = 0 ;
// generazione dei toy experiment e calcolo del fit per ciascuno di essi
// ----------------------------------------
//loop over toys
for (int i_toy = 0 ; i_toy < N_toys ; ++i_toy)
{
vector<double> asse_x ;
vector<double> asse_y ;
// generare il sample
// --------------------
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
double epsilon = rand_TAC_gaus (sigma) ;
asse_x.push_back (i_point) ;
asse_y.push_back (g (i_point) + epsilon) ;
}
// trovare i parametri
// --------------------
matrice H (N_points, 2) ;
for (int i = 0 ; i < N_points ; ++i)
{
H.setCoord (i, 0, 1) ;
H.setCoord (i, 1, asse_x.at (i)) ;
}
vettore y (asse_y) ;
matrice V (N_points) ;
for (int i = 0 ; i < N_points ; ++i) V.setCoord (i, i, sigma * sigma) ;
matrice V_inv = V.inversa () ;
matrice theta_v = (H.trasposta () * V_inv * H).inversa () ;
vettore theta = (theta_v * (H.trasposta () * V_inv)) * y ;
// riempire istogrammi e contatori
// --------------------
h_a.Fill (theta.at (0)) ;
h_b.Fill (theta.at (1)) ;
h_ab.Fill (theta.at (0), theta.at (1)) ;
if (fabs (theta.at (0) - 3.14) < sqrt (theta_v.at (0,0))) ++cont_a ;
if (fabs (theta.at (1) - 2. ) < sqrt (theta_v.at (1,1))) ++cont_b ;
} //loop over toys
cout << "copertura parametro a: " << static_cast<double> (cont_a) / N_toys << endl ;
cout << "copertura parametro b: " << static_cast<double> (cont_b) / N_toys << endl ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetRightMargin (0.15) ;
h_a.SetFillColor (kOrange + 1) ;
h_a.SetLineColor (kGray + 1) ;
h_a.Draw ("hist") ;
c1.Print ("parametro_a.png", "png") ;
h_b.Draw ("hist") ;
h_b.SetFillColor (kOrange + 1) ;
h_b.SetLineColor (kGray + 1) ;
c1.Print ("parametro_b.png", "png") ;
h_ab.Draw ("colz") ;
c1.Print ("parametri_2D.png", "png") ;
return 0 ;
}
A4.3 determinazione della varianza degli scarti delle singole misure y<sub>i</sub>¶
main_03.cpp
/*
c++ -o main_03 `root-config --glibs --cflags` algebra_2.cc main_03.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
#include "TH1F.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TFile.h"
#include "algebra_2.h"
#include "generazione.h"
// retta con termine noto non nullo
double g (double x)
{
return 3.14 + 2 * x ;
}
using namespace std ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv)
{
if (argc < 4)
{
cout << "usage: " << argv[0] << " sigma_y N_points Ntoys" << endl ;
exit (1) ;
}
double sigma = atof (argv[1]) ;
int N_points = atoi (argv[2]) ;
int N_toys = atoi (argv[3]) ;
// istogrammi per il disegno dei risultati del fit
TH1F h_scarti ("h_scarti", "scarti", 200, 0., 5 * N_points) ;
TH1F h_varianza ("h_varianza", "varianza calcolata", 50, 0., 10 * sigma) ;
TH1F h_sigma ("h_sigma", "sigma calcolata", 50, 0., 3 * sigma) ;
// generazione dei toy experiment e calcolo del fit per ciascuno di essi
// ----------------------------------------
//loop over toys
for (int i_toy = 0 ; i_toy < N_toys ; ++i_toy)
{
vector<double> asse_x ;
vector<double> asse_y ;
// generare il sample
// --------------------
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
double epsilon = rand_TAC_gaus (sigma) ;
asse_x.push_back (i_point) ;
asse_y.push_back (g (i_point) + epsilon) ;
}
// trovare i parametri
// --------------------
matrice H (N_points, 2) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
H.setCoord (i_point, 0, 1) ;
H.setCoord (i_point, 1, asse_x.at (i_point)) ;
}
vettore y (asse_y) ;
matrice V (N_points) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
V.setCoord (i_point, i_point, 1.) ;
matrice V_inv = V.inversa () ;
matrice theta_v = (H.trasposta () * V_inv * H).inversa () ;
vettore theta = (theta_v * (H.trasposta () * V_inv)) * y ;
// calcolo la somma degli scarti quadratici normalizzati per l'incertezza
double Q2min = (y - H * theta).dot (V_inv * (y - H * theta)) ;
h_scarti.Fill (Q2min) ;
double varianza = Q2min / (N_points - 2) ;
h_varianza.Fill (varianza) ;
h_sigma.Fill (sigma) ;
} //loop over toys
cout << "sigma = "
<< sqrt (h_scarti.GetMean () / (N_points - 2))
<< endl ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetRightMargin (0.15) ;
h_scarti.Draw ("hist") ;
c1.Print ("scarti.png", "png") ;
TFile f_out ("main_04.root", "recreate") ;
h_scarti.Write () ;
h_sigma.Write () ;
h_varianza.Write () ;
f_out.Close () ;
return 0 ;
}
A4.4 confronto della distribuzione di Q<sup>2</sup><sub>min</sub> con la forma Χ<sup>2</sup>(x, N-k) con scarti non gaussiani¶
main_04.cpp
/*
c++ -o main_04 `root-config --glibs --cflags` algebra_2.cc main_04.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
#include "TH1F.h"
#include "TCanvas.h"
#include "TH2F.h"
#include "TFile.h"
#include "algebra_2.h"
#include "generazione.h"
// retta con termine noto non nullo
double g (double x)
{
return 3.14 + 2 * x ;
}
using namespace std ;
// ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ---- ----
int main (int argc, char ** argv)
{
if (argc < 4)
{
cout << "usage: " << argv[0] << " sigma_y N_points Ntoys" << endl ;
exit (1) ;
}
double sigma = atof (argv[1]) ;
int N_points = atoi (argv[2]) ;
int N_toys = atoi (argv[3]) ;
// istogrammi per il disegno dei risultati del fit
TH1F h_scarti ("h_scarti", "scarti", 200, 0., 5 * N_points) ;
TH1F h_varianza ("h_varianza", "varianza calcolata", 50, 0., 10 * sigma) ;
TH1F h_sigma ("h_sigma", "sigma calcolata", 50, 0., 3 * sigma) ;
// generazione dei toy experiment e calcolo del fit per ciascuno di essi
// ----------------------------------------
//loop over toys
for (int i_toy = 0 ; i_toy < N_toys ; ++i_toy)
{
vector<double> asse_x ;
vector<double> asse_y ;
// generare il sample
// --------------------
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
double epsilon = rand_TAC_uniform (sigma) ;
asse_x.push_back (i_point) ;
asse_y.push_back (g (i_point) + epsilon) ;
}
// trovare i parametri
// --------------------
matrice H (N_points, 2) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
{
H.setCoord (i_point, 0, 1) ;
H.setCoord (i_point, 1, asse_x.at (i_point)) ;
}
vettore y (asse_y) ;
matrice V (N_points) ;
for (int i_point = 0 ; i_point < N_points ; ++i_point)
V.setCoord (i_point, i_point, sigma * sigma) ;
matrice V_inv = V.inversa () ;
matrice theta_v = (H.trasposta () * V_inv * H).inversa () ;
vettore theta = (theta_v * (H.trasposta () * V_inv)) * y ;
// calcolo la somma degli scarti quadratici normalizzati per l'incertezza
double Q2min = (y - H * theta).dot (V_inv * (y - H * theta)) ;
h_scarti.Fill (Q2min) ;
double varianza = Q2min / (N_points - 2) ;
h_varianza.Fill (varianza) ;
h_sigma.Fill (sigma) ;
} //loop over toys
cout << "sigma = "
<< sqrt (h_scarti.GetMean () / (N_points - 2))
<< endl ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetRightMargin (0.15) ;
h_scarti.Draw ("hist") ;
c1.Print ("scarti.png", "png") ;
return 0 ;
}