Esempi: Lezione 12
Contents
Esempi: Lezione 12¶
12.0 Creazione e disegno di un modello funzionale con la classe TF1 di ROOT¶
main_00.cpp
/*
c++ -o main_00 `root-config --glibs --cflags` main_00.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
using namespace std ;
int main (int argc, char ** argv)
{
double N_bkg = 100. ;
double p0 = log (N_bkg) ; // eventi di fondo
double p1 = -0.2 ; // pendenza del fondo
double p2 = 30. ; // eventi di segnale
double p3 = 10. ; // media del segnale
double p4 = 2. ; // sigma del segnale
// The total is the sum of the three, each has 3 parameters
TF1 model ("model", "expo(0) + gaus(2)", 0., 20.) ;
model.SetParameter (0, p0) ;
model.SetParameter (1, p1) ;
model.SetParameter (2, p2) ;
model.SetParameter (3, p3) ;
model.SetParameter (4, p4) ;
model.SetLineColor (kBlue + 2) ;
model.SetLineWidth (4) ;
model.SetLineStyle (1) ;
// The total is the sum of the three, each has 3 parameters
TF1 segnale ("segnale", "gaus(0)", 0., 20.) ;
segnale.SetParameter (0, p2) ;
segnale.SetParameter (1, p3) ;
segnale.SetParameter (2, p4) ;
segnale.SetLineColor (kBlue) ;
segnale.SetLineWidth (2) ;
segnale.SetLineStyle (7) ;
// The total is the sum of the three, each has 3 parameters
TF1 fondo ("fondo", "expo(0)", 0., 20.) ;
fondo.SetParameter (0, p0) ;
fondo.SetParameter (1, p1) ;
fondo.SetLineColor (kGray + 2) ;
fondo.SetLineWidth (2) ;
fondo.SetLineStyle (7) ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
TH1F * bkg = c1.DrawFrame (0, 0, 20, 1.05 * N_bkg) ;
bkg->GetXaxis ()->SetTitle ("x") ;
bkg->GetYaxis ()->SetTitle ("f(x)") ;
segnale.Draw ("same") ;
fondo.Draw ("same") ;
model.Draw ("same") ;
c1.Print ("model.png", "png") ;
// 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 ;
}
12.1 Generazione di numeri pseudo-casuali secondo una distribuzione data tramite ROOT¶
main_01.cpp
/*
c++ -o main_01 `root-config --glibs --cflags` main_01.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 3)
{
cout << "usage: " << argv[0] << " outputfile.txt N_events" << endl ;
exit (1) ;
}
double N_bkg = 100. ;
double p0 = log (N_bkg) ; // eventi di fondo
double p1 = -0.2 ; // pendenza del fondo
double p2 = 30. ; // eventi di segnale
double p3 = 10. ; // media del segnale
double p4 = 2. ; // sigma del segnale
TF1 model ("model", "expo(0) + gaus(2)", 0., 20.) ;
model.SetParameter (0, p0) ;
model.SetParameter (1, p1) ;
model.SetParameter (2, p2) ;
model.SetParameter (3, p3) ;
model.SetParameter (4, p4) ;
TH1F h_campione ("h_campione", "", 100, 0., 20.) ;
ofstream f_campione ;
f_campione.open (argv[1]) ;
int N_eventi = atoi (argv[2]) ;
for (int i = 0 ; i < N_eventi ; ++i)
{
double event = model.GetRandom () ;
f_campione << event << "\n" ;
h_campione.Fill (event) ;
}
f_campione.close () ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
h_campione.GetXaxis ()->SetTitle ("x") ;
h_campione.GetYaxis ()->SetTitle ("eventi nel bin") ;
h_campione.SetFillColor (kOrange + 1) ;
h_campione.SetLineColor (kGray + 1) ;
h_campione.SetStats (0) ;
h_campione.Draw () ;
c1.Print ("sample.png", "png") ;
return 0 ;
}
12.2 Fit di un istogramma costruito con gli eventi generati nell’esercizio precedente¶
main_02.cpp
/*
c++ -o main_02 `root-config --glibs --cflags` main_02.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
using namespace std ;
float minimo (const vector<double> & v)
{
double min = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) < min) min = v.at (i) ;
return min ;
}
float massimo (const vector<double> & v)
{
double max = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) > max) max = v.at (i) ;
return max ;
}
int main (int argc, char ** argv)
{
if (argc < 2)
{
cout << "usage: " << argv[0] << " inputfile.txt" << endl ;
exit (1) ;
}
// lettura del file di eventi
// --------------------------
ifstream f_campione ;
f_campione.open (argv[1]) ;
vector<double> v_eventi ;
while (true)
{
double evento ;
f_campione >> evento ;
if (f_campione.eof () == true) break ;
v_eventi.push_back (evento) ;
}
f_campione.close () ;
cout << "letti " << v_eventi.size () << " eventi" << endl ;
cout << "minimo degli eventi: " << minimo (v_eventi) << endl ;
cout << "massimo degli eventi: " << massimo (v_eventi) << endl ;
// preparazione e riempimento dell'istogramma
// --------------------------
double min = floor (minimo (v_eventi)) ;
double max = ceil (massimo (v_eventi)) ;
int N_bin = round (sqrt (v_eventi.size ()) / 2.) ;
TH1F h_eventi ("h_eventi", "", N_bin, min, max) ;
for (int i = 0 ; i < v_eventi.size () ; ++i) h_eventi.Fill (v_eventi.at (i)) ;
// preparazione del modello per il fit
// --------------------------
TF1 model ("model", "expo(0) + gaus(2)", 0., 20.) ;
model.SetLineColor (kBlue + 2) ;
model.SetLineWidth (4) ;
model.SetLineStyle (1) ;
// determinazione dei parametri iniziali del fit
// --------------------------
// prima stima dei parametri
double N_bkg = v_eventi.size () / 2. ;
double p0 = log (N_bkg) ; // eventi di fondo
double p1 = -0.5 ; // pendenza del fondo
double p2 = v_eventi.size () / 2. ; // eventi di segnale
double p3 = 0.5 * (max - min) ; // media del segnale
double p4 = h_eventi.GetRMS () ; // sigma del segnale
// prima stima di p0 e p1 con un fit in zona di solo fondo
TF1 fondo ("fondo", "expo(0)", 0., 20.) ;
fondo.SetParameter (0, p0) ;
fondo.SetParameter (1, p1) ;
h_eventi.Fit ("fondo", "Q", "", 0., 4.) ;
// prima stima di p2, p3 e p4 con un fit in zona di solo fondo
TF1 segnale ("segnale", "gaus(0)", 0., 20.) ;
segnale.SetParameter (0, p2) ;
segnale.SetParameter (1, p3) ;
segnale.SetParameter (2, p4) ;
h_eventi.Fit ("segnale", "Q", "", 7., 14.) ;
model.SetParameter (0, fondo.GetParameter (0)) ;
model.SetParameter (1, fondo.GetParameter (1)) ;
model.SetParameter (2, segnale.GetParameter (0)) ;
model.SetParameter (3, segnale.GetParameter (1)) ;
model.SetParameter (4, segnale.GetParameter (2)) ;
TFitResultPtr fit_result = h_eventi.Fit ("model", "S") ;
// bonta' del fit
// --------------------------
int result = fit_result ;
cout << "convergenza del fit : " << fit_result->IsValid () << endl ;
cout << "convergenza del fit (bis) : " << fit_result->Status () << endl ;
fit_result->Print () ;
fit_result->PrintCovMatrix (cout) ;
cout << "probabilità associata a Q2 : " << fit_result->Prob () << endl ;
cout << "Valore di Q2 : " << fit_result->Chi2 () << endl ;
cout << "Numero di gradi di libertà : " << fit_result->Ndf () << endl ;
// output dei risultati
// --------------------------
cout << endl ;
cout.precision (3) ;
cout << "normalizzazione del fondo : " << exp (model.GetParameter (0)) << "\t+- "
<< model.GetParError (0) * exp (model.GetParameter (0)) << endl ;
cout << "pendenza del fondo : " << model.GetParameter (1) << "\t+- " << model.GetParError (1) << endl ;
cout << "normalizzazione del segnale: " << model.GetParameter (2) << "\t+- " << model.GetParError (2) << endl ;
cout << "media del segnale : " << model.GetParameter (3) << "\t+- " << model.GetParError (3) << endl ;
cout << "sigma del segnale : " << model.GetParameter (4) << "\t+- " << model.GetParError (4) << endl ;
// matrice di covarianza dei parametri
// --------------------------
cout << endl ;
TMatrixDSym cov = fit_result->GetCovarianceMatrix () ;
// or TMatrixDSym cov = r->GetCorrelationMatrix();
for (int i = 0; i < cov.GetNrows () ; ++i)
{
for (int j = 0; j < cov.GetNcols () ; ++j)
{
cout << cov(i,j) << "\t" ;
}
cout << "\n";
}
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
h_eventi.GetXaxis ()->SetTitle ("x") ;
h_eventi.GetYaxis ()->SetTitle ("eventi nel bin") ;
h_eventi.SetFillColor (kOrange + 1) ;
h_eventi.SetLineColor (kGray + 1) ;
h_eventi.SetStats (0) ;
h_eventi.Draw () ;
c1.Print ("main_02.png", "png") ;
return 0 ;
}
12.3 Generazione di numeri pseudo-casuali secondo una distribuzione Gaussiana tramite ROOT¶
main_03.cpp
/*
c++ -o main_03 `root-config --glibs --cflags` main_03.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
using namespace std ;
int main (int argc, char ** argv)
{
if (argc < 3)
{
cout << "usage: " << argv[0] << " outputfile.txt N_events" << endl ;
exit (1) ;
}
double p0 = 30. ; // integrale
double p1 = 10. ; // media
double p2 = 2. ; // sigma
TF1 * model = new TF1 ("model", "gaus(0)", 0., 20.) ;
model->SetParameter (0, p0) ;
model->SetParameter (1, p1) ;
model->SetParameter (2, p2) ;
TH1F h_campione ("h_campione", "", 100, 0., 20.) ;
ofstream f_campione ;
f_campione.open (argv[1]) ;
int N_eventi = atoi (argv[2]) ;
for (int i = 0 ; i < N_eventi ; ++i)
{
double event = model->GetRandom () ;
f_campione << event << "\n" ;
h_campione.Fill (event) ;
}
f_campione.close () ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
h_campione.GetXaxis ()->SetTitle ("x") ;
h_campione.GetYaxis ()->SetTitle ("eventi nel bin") ;
h_campione.SetFillColor (kOrange + 1) ;
h_campione.SetLineColor (kGray + 1) ;
h_campione.SetStats (0) ;
h_campione.Draw () ;
c1.Print ("gauss.png", "png") ;
delete model ;
return 0 ;
}
12.4 Confronto fra un fit effettuato con i minimi quadrati e con la massima verosimiglianza¶
main_04.cpp
/*
c++ -o main_04 `root-config --glibs --cflags` main_04.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
using namespace std ;
double minimo (const vector<double> & v)
{
double min = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) < min) min = v.at (i) ;
return min ;
}
double massimo (const vector<double> & v)
{
double max = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) > max) max = v.at (i) ;
return max ;
}
double maxval (double a, double b)
{
if (a > b) return a ;
return b ;
}
int main (int argc, char ** argv)
{
if (argc < 2)
{
cout << "usage: " << argv[0] << " inputfile.txt" << endl ;
exit (1) ;
}
// lettura del file di eventi
// --------------------------
ifstream f_campione ;
f_campione.open (argv[1]) ;
vector<double> v_eventi ;
while (true)
{
double evento ;
f_campione >> evento ;
if (f_campione.eof () == true) break ;
v_eventi.push_back (evento) ;
}
f_campione.close () ;
cout << "letti " << v_eventi.size () << " eventi" << endl ;
cout << "minimo degli eventi: " << minimo (v_eventi) << endl ;
cout << "massimo degli eventi: " << massimo (v_eventi) << endl ;
// preparazione e riempimento dell'istogramma
// --------------------------
double min = floor (minimo (v_eventi)) ;
double max = ceil (massimo (v_eventi)) ;
// int N_bin = round (sqrt (v_eventi.size ()) / 2.) ;
int N_bin = round (sqrt (v_eventi.size ())) ;
N_bin = maxval (N_bin, 10) ;
TH1F h_eventi ("h_eventi", "", N_bin, min, max) ;
for (int i = 0 ; i < v_eventi.size () ; ++i) h_eventi.Fill (v_eventi.at (i)) ;
TF1 * model = new TF1 ("model", "gaus(0)", 0., 20.) ;
model->SetLineColor (kBlue + 2) ;
model->SetLineWidth (4) ;
model->SetLineStyle (1) ;
double p0 = v_eventi.size () ; // eventi
double p1 = 0.5 * (max - min) ; // media
double p2 = h_eventi.GetRMS () ; // sigma
model->SetParameter (0, p0) ;
model->SetParameter (1, p1) ;
model->SetParameter (2, p2) ;
TFitResultPtr fit_result_MQ = h_eventi.Fit ("model", "SQ+") ;
h_eventi.GetFunction ("model")->SetLineColor (kRed) ;
cout << endl ;
cout.precision (3) ;
cout << "probabilità associata a Q2: " << model->GetProb () << endl ;
cout << "integrale : " << model->GetParameter (0) << "\t+- " << model->GetParError (0) << endl ;
cout << "media : " << model->GetParameter (1) << "\t+- " << model->GetParError (1) << endl ;
cout << "sigma : " << model->GetParameter (2) << "\t+- " << model->GetParError (2) << endl ;
TFitResultPtr fit_result_ML = h_eventi.Fit ("model", "SLQ+") ;
cout << endl ;
cout << "probabilità associata a Q2: " << model->GetProb () << endl ;
cout << "integrale : " << model->GetParameter (0) << "\t+- " << model->GetParError (0) << endl ;
cout << "media : " << model->GetParameter (1) << "\t+- " << model->GetParError (1) << endl ;
cout << "sigma : " << model->GetParameter (2) << "\t+- " << model->GetParError (2) << endl ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
h_eventi.GetXaxis ()->SetTitle ("x") ;
h_eventi.GetYaxis ()->SetTitle ("eventi nel bin") ;
h_eventi.SetFillColor (kOrange + 1) ;
h_eventi.SetLineColor (kGray + 1) ;
h_eventi.SetStats (0) ;
h_eventi.Draw () ;
c1.Print ("main_04.png", "png") ;
return 0 ;
}
12.5 Studio del confronto in funzione del numero di eventi¶
main_05.cpp
/*
c++ -o main_05 `root-config --glibs --cflags` main_05.cpp
*/
#include <cstdlib>
#include <iostream>
#include <cmath>
#include <fstream>
#include "TF1.h"
#include "TH1F.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TFitResult.h"
#include "TMatrixDSym.h"
#include "TGraphErrors.h"
using namespace std ;
double minimo (const vector<double> & v)
{
double min = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) < min) min = v.at (i) ;
return min ;
}
double massimo (const vector<double> & v)
{
double max = v.at (0) ;
for (int i = 0 ; i < v.size () ; ++i)
if (v.at (i) > max) max = v.at (i) ;
return max ;
}
double maxval (double a, double b)
{
if (a > b) return a ;
return b ;
}
template <class disegno>
void drawComparison (TCanvas & c1, disegno & d1, disegno & d2, TString variabile)
{
d1.SetMarkerStyle (20) ;
d2.SetMarkerStyle (4) ;
d1.SetMarkerColor (kBlue) ;
d2.SetMarkerColor (kRed) ;
d1.Draw ("AP") ;
d1.GetXaxis ()->SetTitle ("numero di eventi") ;
d1.GetYaxis ()->SetTitle (variabile) ;
d1.Draw ("AP") ;
d2.Draw ("P same") ;
c1.Print (variabile + ".png", "png") ;
return ;
}
int main (int argc, char ** argv)
{
if (argc < 2)
{
cout << "usage: " << argv[0] << " inputfile.txt" << endl ;
exit (1) ;
}
// lettura del file di eventi
// --------------------------
ifstream f_campione ;
f_campione.open (argv[1]) ;
vector<double> v_eventi ;
while (true)
{
double evento ;
f_campione >> evento ;
if (f_campione.eof () == true) break ;
v_eventi.push_back (evento) ;
}
f_campione.close () ;
// preparazione dei TGraph per l'andamento dei risultati
// --------------------------
vector<TGraphErrors *> trend_parametri ;
for (int i = 0 ; i < 5 ; ++i) trend_parametri.push_back (new TGraphErrors ()) ;
vector<TGraph *> trend_incertezza_par ;
for (int i = 0 ; i < 5 ; ++i) trend_incertezza_par.push_back (new TGraph ()) ;
TGraph g_gof ;
TCanvas c1 ("c1", "", 800, 800) ;
c1.SetLeftMargin (0.15) ;
TGraphErrors g_p0_MQ ;
TGraphErrors g_p1_MQ ;
TGraphErrors g_p2_MQ ;
TGraphErrors g_p0_ML ;
TGraphErrors g_p1_ML ;
TGraphErrors g_p2_ML ;
int point = 0 ;
// ciclo su diverse dimensioni del campione
for (int i_Nmax = v_eventi.size () ; i_Nmax > 10 ; i_Nmax /= 2)
{
double min = floor (minimo (v_eventi)) ;
double max = ceil (massimo (v_eventi)) ;
int N_bin = round (sqrt (i_Nmax)) ;
N_bin = maxval (N_bin, 10) ;
TH1F * h_eventi = new TH1F ("h_eventi", "", N_bin, min, max) ;
for (int i = 0 ; i < i_Nmax ; ++i) h_eventi->Fill (v_eventi.at (i)) ;
// preparazione del modello per il fit
// --------------------------
TF1 * model = new TF1 ("model", "gaus(0)", 0., 20.) ;
model->SetLineColor (kBlue + 2) ;
model->SetLineWidth (4) ;
model->SetLineStyle (1) ;
double p0 = i_Nmax ; // eventi
double p1 = 0.5 * (max - min) ; // media
double p2 = h_eventi->GetRMS () / 2.; // sigma
model->SetParameter (0, p0) ;
model->SetParameter (1, p1) ;
model->SetParameter (2, p2) ;
TFitResultPtr fit_result_MQ = h_eventi->Fit ("model", "SQ+") ;
h_eventi->GetFunction ("model")->SetLineColor (kRed) ;
g_p0_MQ.SetPoint (point, i_Nmax, model->GetParameter (0)) ;
g_p0_MQ.SetPointError (point, 0., model->GetParError (0)) ;
g_p1_MQ.SetPoint (point, i_Nmax, model->GetParameter (1)) ;
g_p1_MQ.SetPointError (point, 0., model->GetParError (1)) ;
g_p2_MQ.SetPoint (point, i_Nmax, model->GetParameter (2)) ;
g_p2_MQ.SetPointError (point, 0., model->GetParError (2)) ;
TFitResultPtr fit_result_ML = h_eventi->Fit ("model", "SLQ+") ;
g_p0_ML.SetPoint (point, i_Nmax, model->GetParameter (0)) ;
g_p0_ML.SetPointError (point, 0., model->GetParError (0)) ;
g_p1_ML.SetPoint (point, i_Nmax, model->GetParameter (1)) ;
g_p1_ML.SetPointError (point, 0., model->GetParError (1)) ;
g_p2_ML.SetPoint (point, i_Nmax, model->GetParameter (2)) ;
g_p2_ML.SetPointError (point, 0., model->GetParError (2)) ;
h_eventi->GetXaxis ()->SetTitle ("x") ;
h_eventi->GetYaxis ()->SetTitle ("eventi nel bin") ;
h_eventi->SetFillColor (kOrange + 1) ;
h_eventi->SetLineColor (kGray + 1) ;
TString testo = "fit_" ;
testo += point ;
testo += ".png" ;
h_eventi->Draw () ;
c1.Print (testo, "png") ;
delete h_eventi ;
delete model ;
++point ;
} // ciclo su diverse dimensioni del campione
c1.SetLogx () ;
drawComparison (c1, g_p0_MQ , g_p0_ML , "p0") ;
drawComparison (c1, g_p1_MQ , g_p1_ML , "p1") ;
drawComparison (c1, g_p2_MQ , g_p2_ML , "p2") ;
return 0 ;
}