Lezione 4
In questa lezione proveremo ad applicare gli strumenti sviluppati nelle scorse lezioni a casi concreti di studio. Per prima cosa cercheremo di chiudere il capitolo dell'analisi dei dati climatici studiando l'andamento temporale delle variazioni medie di temperatura dal 1941 al 2023. Vi proponiamo inoltre alcuni temi specifici (facoltativi) che si possono affrontare con gli strumenti che abbiamo sviluppato sinora: l'analisi dei dati raccolti in un tipico esperimento del laboratorio di fisica (esperimento di Millikan e misura del rapporto q/m per l'elettrone) e un semplice problema di riordinamento (per fare qualche esercizio con gli strumenti della STL).
ESERCIZIO 4.0 - Analisi andamento temporale temperature Milano (da consegnare)
In questa lezione tireremo le somme del lavoro svolto in tutte le sessioni precedenti. Scriveremo un codice di analisi (da consegnare) che ci permetta di visualizzare l'andamento del delta giornaliero (calcolato rispetto alla media sul perioro 1941-2023) medio annuale rispetto all'anno dal 1941 al 2023 nell'area di Milano. Per analizzare l'andamento temporale del in funzione del tempo l'oggetto giusto di ROOT da utilzzare è il TGraphErrors. La struttura del codice potrebbe essere logicamente la seguente :
-
Loop principale sugli anni che vogliamo investigare
-
Per ogni anno apriamo il file corrispondente, calcoliamo media ed errore e riempiamo un grafico
-
A fine ciclo disegnamo il grafico
I files di dati li potete scaricare qui.
Esempio di codice :
// put here all required includes
using namespace std;
int main( ) {
TApplication app("app",0,0);
// oggetto per rappresentare un andamento x (anno) verso y (delta medio)
TGraphErrors trend;
// ...
int index {0} ;
// loop principale su tutti i files ( anni ) da analizzare
for ( int i = 1941 ; i < 2024 ; i++) {
string filename = to_string(i) + ".txt" ;
vector<double> v = Read<double> ( filename.c_str() ) ;
// qui fate i vostri calcoli
cout << " Anno " << to_string(i) << " delta medio = " << ave << " +/- " << err << endl;
// inserisco media e deviazione standard dalla media nel grafico
trend.SetPoint(index, i, ave);
trend.SetPointError( index , 0 , err);
index++;
}
// grafica
TCanvas c("Temperature trend","Temperature trend");
c.cd();
c.SetGridx();
c.SetGridy();
trend.SetMarkerSize(1.0);
trend.SetMarkerStyle(20);
trend.SetFillColor(5);
trend.SetTitle("Temperature trend");
trend.GetXaxis()->SetTitle("Anno");
trend.GetYaxis()->SetTitle("#Delta (#circ C)");
trend.Draw("apl3");
trend.Draw("pX");
c.SaveAs("trend.pdf");
app.Run();
}
Quale errore associare ai ?
In linea di principio lo stimatore corretto sarebbe la deviazione standard dalla media. In questo caso la stima è più complicata perchè i valori di giornalieri non sono scorrelati tra loro. Per ridurre l'impatto della correlazione e quindi utilizzare uno stimatore dell'errore più corretto potremmo limitaci ad utilizzare una misura ogni 7 giorni: ci aspettiamo infatti che una misura sia solo debolmente correlata a quello che è successo 7 giorni prima.
Per svolgere questo esercizio (e tutti i successivi) potete chiaramente utilizzare lo strumento di rappresentazione grafica che preferite (GNUPLOT o Matplotlib).
ESERCIZIO 4.1 - Misura del rapporto q/m per l'elettrone
Un esperimento molto interessante che si svolge nel laboratorio di fisica riguarda la misura del rapporto tra la carica e la massa di un elettrone. Un fascio di elettroni di energia nota (determinata dal potenziale 2ΔV applicato ai capi di un condensatore) viene deflesso da un campo magnetico opportunamente generato. Dalla misura del raggio di deflessione r si può ricavare una stima del rapporto q/m. La misura in laboratorio è piuttosto complessa ( in particolare lo studio delle sistematiche sperimentali ) ma il rapporto q/m può essere determinato sfruttando la relazione
come il coefficiente angolare della relazione. Nel file di dati sono riportate le misure prese in laboratorio: nella prima colonna il valore di , nella seconda colonna il valore di e nella terza l'errore sulla determinazione di . Provate a scrivere un codice che pemetta di visualizzare i dati raccolti e determinare il rapporto q/m ( in questo caso l'oggetto di ROOT più indicato è il TGraphErrors ).
- Il file di dati lo potete trovare qui.
- Potete provare a scrivere il codice da soli. In caso potete prendere ispirazione dall'esempio qui sotto.
Esempio di codice :
#include <iostream>
#include <vector>
#include <fstream>
#include "TGraphErrors.h"
#include "TApplication.h"
#include "TF1.h"
#include "TAxis.h"
#include "TLegend.h"
#include <cmath>
using namespace std;
void ParseFile( string filename, vector<double> & myx, vector<double> & myy , vector<double> & myerry ) {
ifstream fin(filename.c_str());
double x, y, err ;
if ( !fin ) {
cout << "Cannot open file " << filename << endl;
exit(11);
};
while ( fin >> x >> y >> err ) {
myx.push_back(x);
myy.push_back(y);
myerry.push_back(err);
};
fin.close();
}
TGraphErrors DoPlot( vector<double> myx, vector<double> myy , vector<double> myerry ) {
TGraphErrors mygraph;
for ( int k = 0 ; k < myx.size() ; k++ ) {
mygraph.SetPoint(k, myx[k] , myy[k] );
mygraph.SetPointError(k, 0, myerry[k]);
};
return mygraph;
}
int main() {
TApplication app("MyApp",0,0);
vector<double> myx ;
vector<double> myy ;
vector<double> myerry;
// read data from file
ParseFile("data_eom.dat", myx, myy, myerry ) ;
// create TGraphErrors
TGraphErrors mygraph = DoPlot(myx, myy, myerry);
// fit the TGraphErrors ( linear fit )
TF1 * myfun = new TF1 ("fitfun","[0]*x+[1]",0, 1000) ;
mygraph.Fit(myfun);
double moe = myfun->GetParameter(0);
double error = sqrt( pow (1./moe, 4) * pow(myfun->GetParError(0),2) ) ;
cout << "Valore di e/m dal fit = " << 1./moe << "+/-" << error << endl;
cout << "Valore del chi2 del fit = " << myfun->GetChisquare() << endl;
cout << " prob del chi2 = " << myfun->GetProb() << endl;
// customise the plot, cosmetics
mygraph.Draw("AP");
mygraph.SetMarkerStyle(20);
mygraph.SetMarkerSize(1);
mygraph.SetTitle("Misura e/m");
mygraph.GetXaxis()->SetTitle("2#DeltaV (V)");
mygraph.GetYaxis()->SetTitle("(B_{z}R)^{2} (T^{2}m^{2} )");
TLegend leg (0.15,0.7, 0.3, 0.85) ;
leg.AddEntry(&mygraph,"data","p");
leg.AddEntry(myfun,"fit","l");
leg.Draw("same");
app.Run();
return 0;
}
ESERCIZIO 4.2 - Misura della carica dell'elettrone
La carica dell'elettrone è stata misurata per la prima volta nel 1909 in uno storico esperimento dal fisico statunitense Robert Millikan. Questo esperimento si effettua anche nel laboratorio di fisica e consiste nel misurare la velocità di caduta o risalita di goccioline d'olio elettrizzate per strofinio in una regione con campo elettrico regolabile. Nel file di dati sono riportare le misure di carica elettrica (Qi) per un certo numero di goccioline osservate. Il valore della carica può essere determinato come il minimo della funzione dove è il numero intero più vicino al rapporto . Provate a scrivere un codice per rappresentare la funzione e determinare il valore della carica dell'elettrone.
- Il file di dati lo potete trovare qui.
- Potete provare a scrivere il codice da soli. In caso potete prendere ispirazione dall'esempio qui sotto.
Esempio di codice :
#include <iostream>
#include <vector>
#include <cmath>
#include <fstream>
#include <iomanip>
#include "TH1F.h"
#include "TApplication.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TAxis.h"
using namespace std;
// ===========================================================
// Read data from file and store them in a vector
// ===========================================================
vector<double> ParseFile( string filename ) {
ifstream fin(filename.c_str());
vector<double> v;
double val;
if ( !fin ) {
cout << "Cannot open file " << filename << endl;
exit(11);
};
while ( fin >> val ) v.push_back(val);
fin.close();
return v;
}
// ===========================================================
// Compute S(q)
// ===========================================================
double fun ( double q , vector<double> params ) {
double sum = 0;
for ( int k = 0 ; k < params.size() ; k++ ) sum+= pow ( q - params[k]/(round(params[k]/q)) , 2 ) ;
return sum;
}
// ===========================================================
// Compute qmin
// ===========================================================
double deriv ( double qmin , vector<double> params) {
double sum = 0;
for ( int k = 0 ; k < params.size() ; k++ ) sum+= ( params[k]/ round(params[k]/qmin) );
return sum/params.size() ;
}
// ===========================================================
// This can be used to actually minimse S(q)
// ===========================================================
/*
double funROOT ( double* q , double * params ) {
double sum = 0;
for ( int k = 0 ; k < 64 ; k++ ) sum+= pow ( q[0] - params[k] / (round(params[k]/q[0])), 2 ) ;
return sum;
}
*/
// ===========================================================
// This code estimates the best value of qe from a set of
// measurements (drop charges)
// ===========================================================
int main() {
TApplication app(0,0,0);
// read charges from file
vector<double> charges = ParseFile("data_millikan.dat");
// show charges distribution
TCanvas can1 ;
can1.cd();
TH1F histo("cariche","Charges distribution", 100, 0 , 20E-19);
for ( auto i = 0; i < charges.size() ; i++ ) histo.Fill(charges[i]);
histo.Draw();
histo.GetXaxis()->SetTitle("Charge [C]");
TGraph g ;
int counter = 0;
double qmin = 0;
double sqmin = DBL_MAX;
for ( double value = 1.4e-19 ; value < 1.8E-19 ; value+=0.001E-19 ) {
g.SetPoint(counter, value,fun( value , charges ) );
if ( fun( value , charges ) < sqmin ) { sqmin = fun( value , charges ) ; qmin = value ; } ;
counter++;
}
cout << "Found approximate minimum at q = " << qmin << endl;
TCanvas can2;
can2.cd();
g.Draw("ALP");
g.SetMarkerStyle(20);
g.SetMarkerSize(0.5);
g.SetTitle("Best charge value");
g.GetXaxis()->SetTitle("Charge (C)") ;
g.GetYaxis()->SetTitle("S(q) (C^{2})") ;
/*
TF1 myfun( "charge",funROOT,1.5e-19 , 1.7E-19, 64);
for ( int k = 0 ; k < 64; k++ ) myfun.SetParameter(k, charges[k]);
cout << myfun.GetMinimumX() << endl ;
*/
double mycharge = deriv(qmin, charges ) ;
double uncer = sqrt( fun(mycharge, charges) / ( charges.size() * (charges.size()-1)) );
cout << "Measured charge = " << mycharge << " +/- " << uncer << "(stat only)" << endl;
app.Run();
}
ESERCIZIO 4.3 - Determinazione del cammino minimo
In questo esercizio possiamo provare ad approfondire l'uso di contenitori e algoritmi della STL. Proviamo ad affrontare questo problema: dato un set di punti nel piano cerchiamo il percorso che ci permette ( partendo dall'origine ) di toccare tutti i punti percorrendo la minor distanza possibile. Le coordinate dei punti si trovano in un file. Eventualmente cercate di produrre un grafico del percorso effettuato.
- Il file di dati lo potete trovare qui.
- Potete provare a scrivere il codice da soli. In caso potete prendere ispirazione dall'esempio qui sotto che implementa un algoritmo in stile brute force (nel codice si può anche trovare un esempio di overloading di operator<< utile per semplificare la lettura da file ).
Esempio di codice :
#include <vector>
#include <iostream>
#include <fstream>
#include <cmath>
#include <algorithm>
#include "TApplication.h"
#include "TGraph.h"
#include "TAxis.h"
#include <string>
using namespace std;
// implementazione di una classe posizione
class Posizione {
public:
Posizione() { m_x = 0 ; m_y = 0 ; m_z = 0; } ;
Posizione( double x, double y, double z ) { m_x = x; m_y=y ; m_z=z ; } ;
friend std::istream& operator>>(std::istream& is , Posizione& p ) {
string temp;
std::getline(is, temp, ',' );
p.m_x = std::stod(temp);
std::getline(is, temp, ',' );
p.m_y = std::stod(temp);
std::getline(is, temp, '\n' );
p.m_z = std::stod(temp);
return is;
};
double getDistance() const { return sqrt( m_x*m_x + m_y * m_y + m_z * m_z ) ; } ;
double GetX() const { return m_x ; } ;
double GetY() const { return m_y ; } ;
double GetZ() const { return m_z ; } ;
double getDistance ( Posizione p ) const {
double dx = p.GetX() - m_x ;
double dy = p.GetY() - m_y ;
double dz = p.GetZ() - m_z ;
return sqrt( dx*dx + dy*dy + dz*dz ) ;
} ;
bool operator< ( const Posizione & b ) const { return ( getDistance() > b.getDistance() ); } ;
void printPositions() { cout << "Posizione : x = " << m_x << " y = " << m_y << " z = " << m_z << endl; };
private:
double m_x, m_y, m_z ;
};
// algoritmo di riordinamento dei vettore di punti
template <typename T> void findBestPath( T start , T end ) {
Posizione ref ( 0, 0,0 );
for ( auto it = start ; it != end ; it++ ) {
sort( it , end , [&] (Posizione i , Posizione j ) { return i.getDistance(ref) < j.getDistance(ref) ; } ) ;
ref = *it ;
}
}
// main
int main() {
TApplication app("myapp",0,0);
ifstream f ( "data_points.dat") ;
vector<Posizione> vp;
Posizione p ;
while ( f.good() ) {
f>>p ;
vp.push_back( p );
}
for ( auto it = vp.begin() ; it != vp.end() ; it++ ) it->printPositions() ;
findBestPath<vector<Posizione>::iterator>( vp.begin() , vp.end() ) ;
for ( auto it = vp.begin() ; it != vp.end() ; it++ ) it->printPositions() ;
TGraph mygraph;
mygraph.SetPoint(0, 0,0 );
int counter = 1;
for ( auto it = vp.begin() ; it != vp.end() ; it++ ) {
mygraph.SetPoint(counter, (*it).GetX(), (*it).GetY() );
counter++;
}
mygraph.GetXaxis()->SetLimits(-10,10);
mygraph.SetMinimum(-10);
mygraph.SetMaximum(10);
mygraph.GetXaxis()->SetTitle("X");
mygraph.GetYaxis()->SetTitle("Y");
mygraph.SetTitle("Percorso");
mygraph.SetMarkerStyle(21);
mygraph.SetMarkerSize(1);
mygraph.SetLineColor(6);
mygraph.Draw("ALP");
app.Run();
}