Laboratorio di Calcolo 2

proff. A. Andreazza, D. Galli, E. Spoletini

Universita' degli Studi di Milano

Anno Accademico 2006/2007

Lezione 6

Introduzione

Scopo di questa sessione di laboratorio e' sviluppare la simulazione di uno degli esperimenti del laboratorio del secondo anno.

Lo sviluppo di simulazioni trova applicazioni in molti problemi di fisica, sia sperimentale che teorica.

Ad esso sara` abbinato l'approfondimento dell'utilizzo dei pacchetti di ROOT per la determinazione dei parametri di una funzione dal best-fit ai dati presenti in un istogramma.

La simulazione di un esperimento

Una misura sperimentale consiste nella raccolta di un insieme di dati, che non sono esatti, ma di cui è generalmente nota la distribuzione statistica (tipicamente si assume una distribuzione gaussiana attorno al valore vero). Generalmente questi dati vengono poi elaborati per dedurne il valore di grandezze fisiche interessanti. Per esempio, dalle misure del periodo di oscillazione $T$ e della lunghezza $l$ di un pendolo, si può ricavare il valore dell'accelerazione di gravità:

$\begin{displaymath} g=\frac{4\pi ^{2}l}{T^{2}}.\end{displaymath}$

In un caso come questo, in cui la dipendenza funzionale della grandezza fisica significativa dalle quantità misurate è relativamente semplice, possiamo stimare l'errore su $g$ attraverso le semplici formule di propagazione degli errori:

$\begin{displaymath} \left( \frac{\sigma _{g}}{g}\right) ^{2}=\left( \frac{\sigma _{l}}{l}\right) ^{2}+\left( 2\frac{\sigma _{T}}{T}\right) ^{2}.\end{displaymath}$

Tuttavia molto spesso un esperimento richiede un'analisi dei dati molto più complessa e non è possibile dare una semplice formula matematica per connettere le variabili misurate (e le loro incertezze) con le variabili derivate (e le loro incertezze). In altri casi l'ipotesi stessa di errori gaussiani non è valida. In tal caso si deve procedere ad una simulazione del processo di analisi dei dati per essere in grado di valutare se il procedimento è corretto e dare una stima delle incertezze.

In pratica il processo logico che si segue e` il seguente:

si assume un valore vero delle costanti fisiche di interesse;
si calcola come questo valore vero si traduce un valore vero delle grandezze osservabili;
si passa dal valore vero delle osservabili a dei possibili valori misurati di tali quantita`, aggiungendo degli errori di misura casuali;
si applica ai valori misurati la stessa procedura che si applicherebbe alle misure fatte per ottenere dei valori misurati delle costanti fisiche che si cercano.

Le domande a cui si vuole rispondere sono:

i valori misurati cosi` calcolati sono una valutazione corretta dei valori veri delle costanti fisiche?
qual e` l'incertezza di questi valori misurati?

Tale processo viene anche utilizzato durante il progetto di nuovi esperimenti, per rispondere a domande del tipo: "Se voglio ottenere un certo errore sulla grandezza fondamentale X, con quale errore devo misurare le quantità osservabili W, Y e Z?"

Tutte queste domande sono di carattere statistico e, mentre spesso l'esperimento possiamo eseguirlo una volta sola, una simulazione puo` essere ripetuta un numero grande di volte e quindi ricavare le proprieta` statistiche dalla distribuzione di risultati.

Errore ed incertezza

In gergo spesso le parole errore ed incertezza sono utilizzate in maniera interscambiabile.

In quanto segue, useremo i due termini in maniera molto piu` specifica (notate che questa non e` una convenzione universale):

il termine errore indica la differenza tra un valore misurato ed il valore vero: se effettuiamo tante misure di una stessa quantita` con uno stesso strumento, otterremo errori diversi ogni volta;
il termine incertezza indica lo scarto quadratico medio della distribuzione degli errori, questa e` una proprieta` dell'apparato sperimentale, per cui se ripetiamo una stessa misura molte volte, l'incertezza di tutte le misure e' identica.

In una simulazione noi sappiamo direttamente l'errore, e dalla sua disctibuzione calcoliamo l'incertezza della misura reale.

In una misura reale, in cui il valore vero non e' conosciuto, possiamo solo parlare di incertezza.

L'esperimento dello spettrometro a prisma

La dipendenza dell'indice di rifrazione dalla lunghezza d'onda della luce incidente viene descritta dalla legge di Cauchy:

L'esperienza dello spettrometro a prisma si propone di misurare l'indice di rifrazione del materiale di un prisma per le diverse lunghezze d'onda di una lampada al mercurio onde determinare i parametri A e B che caratterizzano tale materiale.

Nell'esperimento l'unico tipo di grandezze che vengono misurate sono gli angoli, per cui possiamo assumere un'incertezza uguale per tutte le misure angolari e pari a:

σ_θ=1' ovvero 0.3 mrad.

Come valori veri, prenderemo invece:

angolo di apertura del prisma α= 50^o
A=2.7
B=60000 nm²

Oppure, chi ha gia` effettuato l'esperimento in laboratorio, puo` utilizzare i valori misurati durante l'esperienza.

L'esperienza di laboratorio prevede la sequenza di misure:

Determinazione dell'angolo di apertura del prisma
La prima misura da fare in laboratorio e' la determinazione dell'angolo di apertura del prisma. Questo si effettuta utilizzando il fascio riflesso.
Supponiamo di aver trovato una posizione in cui vediamo il fascio riflesso da una faccia del prisma. Se a questo punto fissiamo lampada e oculare, possiamo poi ruotare il prisma in modo tale che la seconda faccia si trovi nello stesso posto della prima e la luce venga riflessa nella stessa direzione.

Figure 1: Misura dell'angolo di apertura del prisma attraverso la riflessione.

$\resizebox*{1\columnwidth}{!}{\includegraphics{prisma.eps}}$

Per portarci in questa posizione, il prisma deve venire ruotato di un angolo pari a 180^o-α, quindi se θ₁ e θ₂ indicano la posizione del prisma rispettivamente prima e dopo la rotazione, abbiamo la relazione:

$\begin{displaymath} \alpha =180^{\circ }-\left( \theta _{1}-\theta _{2}\right) . \end{displaymath}$
Misura del fascio indeflesso
Rimuovendo il prisma dal goniometro, si puo` poi misurare la posizione angolare in cui si osserva il fascio della lampada al mercurio.
Indichiamo con θ₀ la posizione del fascio indeflesso.
Misura degli angoli di deviazione minima
Rimettendo il prisma, si va poi a cercare la posizione del fascio rifratto dal prisma stesso. Una volta individuato il fascio corrispondente ad una certa lunghezza d'onda, proviamo a muovere lentamente il prisma. Il fascio rifratto si spostera` di conseguenza.
Tuttavia, se cerchiamo di spostare il prisma in modo che il fascio rifratto si muova verso la posizione θ₀, vediamo che il fascio raggiunge una distanza minima da questa posizione e poi torna indietro.
Misuriamo la posizione angolare del fascio in questo punto θ_m(λ) per le diverse lunghezze d'onda della lampada al mercurio:

Table 1: Linee nel visibile del tubo di Plucker a vapori di Hg

Linea $\lambda$ [nm]

Giallo 1 579.1

Giallo 2 577.0

Verde 546.1

Azzurro 491.6

Indaco 435.8

Viola 1 407.7

Viola 2 404.7

L'angolo di deviazione minima e` dato da:
$\begin{displaymath} \delta _{m}\left( \lambda \right) =\theta _{m}\left( \lambda \right) -\theta _{0}\end{displaymath}$
Determinazione degli indici di rifrazione
Una volta determinato l'angolo di deflessione minima, questo puo` essere collegato all'indice di rifrazione dalla relazione:

$\begin{displaymath} n\left( \lambda \right) =\frac{\sin \frac{\delta _{m}\left( \lambda \right) +\alpha }{2}}{\sin \frac{\alpha }{2}}. \end{displaymath}$

In questo modo possiamo determinare gli indici di rifrazione n(λ) per le diverse lunghezze d'onda della lampada.
Determinazione dei parametri A e B
Una volta determinati tutti gli indici di rifrazione, si pone il problema di determinare i valori dei parametri A e B della formula di Caucy.
Per analizzare i dati rapidamente in laboratorio si utilizza una linearizzazione della formula, graficando n²(λ) in funzione di 1/λ², ed applicando una regressione lineare:

$\begin{displaymath} n\left( \lambda \right) ^{2}=A+\frac{B}{\lambda ^{2}}. \end{displaymath}$ (1)

Tuttavia, non e` sempre possibile ricavare una relazione lineare in funzione dei parametri che ci interessano. Quindi in laboratorio di calcolo impareremo come effettuare una determinazione di tali parametri usando la formula non lineare

**Figure 1:** Misura dell'angolo di apertura del prisma attraverso la riflessione.
$\resizebox*{1\columnwidth}{!}{\includegraphics{prisma.eps}}$

Il programma di simulazione

Il programma prevede di simulare un esperimento completo:

Determinare i valori veri di θ₁, θ₂, θ₀ e dei vari θ_m(λ)
Aggiungere ai valori veri i relativi errori sperimentali, usando una distribuzione gaussiana con larghezza σ_θ.
Calcolare i valori misurato degli indici di rifrazione n(λ).
Effettuare un fit per la determinazione di A e B.

Questo esperimento deve essere ripetuto 10000 volte, ed i risultati di ogni esperimento stampati su di un file in modo da poter poi essere analizzati con ROOT.

Si noti che la posizione iniziale della lampada e del prisma sono arbitrarie: dipendono da come li ha lasciati il gruppo che ha fatto l'esperienza il giorno prima, o potrebbero essere stati spostati casualmente. Quindi i valori veri di θ₀ e θ₁ sono casuali e diversi da esperimento a esperimento: possono essere presi da una distribuzione uniforme tra 0^o e 360^o.

La struttura del programma potrebbe essere come segue:


int main() {
  int Nesperimenti=10000;
  int Nlambda=7;
  double n=new double[Nlambda];
  // inizializzazione delle altre variabili
  ...
  ofstream f("risultati.dat");
  TRandom generatore(0);
  for ( int i=0; i<Nesperimenti; i++) {
    // calcolo dei valori veri per questo esperimento
    ...
    // aggiunta degli errori di misura
    ...
    for (int j=0; j<Nlambda; j++) {
      n[j]=...  // calcolo dell'indice di rifrazione
      f << "\t" << n[j];
    }
    // per fare il fit che determina A e B, utilizzeremo una funzione
    // particolare, in un primo momento possiamo tralasciare 
    // questo punto che completeremo piu` tardi.
    // parametriCauchy(A,sigmaA,B,sigmaB,n);     
    // f << "\t" << A << "\t" << B;
    f << endl;
  }
  f.close();
  return 0;
}

Ci preoccuperemo in un secondo tempo di realizzare la funzione per determinare A e B dagli indici di rifrazione misurati, per il momento limitiamoci a calcolare gli indici di rifrazione misurati ed osservare delle quantita` relative ad essi.

Compilazione con classi di ROOT

Questa parte non e` fattibile con la distributione di ROOT per Windows distribuita nel CD, quindi deve essere eseguita o collegandosi ai PC del laboratorio, oppure con delle macro, invece che con un programma compilato.

Leggendo il frammento di codice precedente, vi sarete resi conto che questo e' un programma compilato (utilizza una funzione main), all'interno del quale vogliamo utilizzare delle classi di ROOT (il generatore di numeri casuali TRandom).

Il C++ e' molto potente nel permetterci di riutilizzare oggetti scritti da altri, ma per poterlo fare questo dobbiamo seguire strettamente alcune regole:

Includere nel codice sorgente l'header file che contiene la definizione degli oggetti che vi interessano:
#include "TRandom.h"
al momento della creazione del file oggetto, specificare al compilatore in quale directory deve andare a cercare gli header file che gli servono. Questo si fa utilizzando l'opzione -Idirectory del g++. Questo si puo` ottenere ad esempio scrivendo nel makefile:
```
%.o : %.cxx
       g++ -c -I$ROOT_SYS/include $<
```
al momento della creazione dell'eseguibile (ovvero al link), bisogna specificare di utilizzare le librerie di ROOT. Una serie di librerie aggiuntive si fornisce al compilatore utilizzando la coppia di opzioni -L e -l:
-Ldirectory -llibreria1 -llibreria2...
Purtroppo la lista delle librerie di ROOT e` piuttosto lunga. Fortunatamente, ROOT ci fornisce un comando che permette di fare automaticamente tale lista:
root-config --libs
Nel makefile, possiamo per comodita` definire una variabile ausiliaria che contiene output del comando:
LIBS= `root-config --libs`
(se includiamo un comando tra due caratteri di accento grave, la shell sostituira` al comando una stringa contenente il suo risultato), e poi usare questa variabile nelle istruzioni di compilazione. Ad esempio, se il nostro int main() fosse in un file simulazione.cxx, avremmo nel makefile la regola:
```
simulazione : simulazione.o
       g++ -Wall -o $@ $^ $LIBS
```

Proprieta` statistiche degli errori sugli n

Per le diverse lunghezze d'onda, fare degli istogrammi degli errori sugli n(λ):

calcolare l'incertezza della misura (se possibile, confrontarla con quella che viene dalla propagazione degli errori)
verificare che il valore medio della distribuzione sia compatibile con 0 (N.B.: l'incertezza sul valor medio non e' l'RMS dell'istogramma ma e` collegata ad essa da una relazione statistica che abbiamo appena ripassato...)

Fare un grafico bidimensionale degli errori per una lunghezza d'onda in funzione degli errori su di un'altra lunghezza d'onda di vostra scelta. Questo grafico dovrebbe mostrare un certo grado di correlazione:

Avere dati correlati, significa in pratica che, se uno dei due e` sbagliato in una certa direzione (troppo alto o troppo basso), anche l'altro tendera` ad avere un errore nella stessa direzione. La correlazione viene valutata quantitativamente calcolando il coefficiente di correlazione:

Calcolatene il valore in questo caso.

Per verificare l'origine della correlazione, provate ad azzerare tutti gli errori e poi a riattivare a turno una sola sorgente di errore. Per ognuna delle quantita` misurate, θ₁, θ₂, θ₀ e θ_m(λ), valutare:

che contributo da` all'incertezza della misura,
se e' una sorgente di correlazione o meno.

Determinazione di A e B

Avendo misurato dei valori di n(λ) per diverse lunghezze d'onda, possiamo definire come migliore stima dei parametri A e B, date le nostre misure sperimentali, quella coppia di parametri che minimizza il χ²:

Ci sono diversi pacchetti di software che effettuano questo tipo di minimizzazioni, noi vedremo come farlo dentro ROOT.

Passo 1: creare un grafico dei dati

Il modo standard e' quello di mettere i dati in un grafico, e poi utilizzare il metodo Fit della classe TGraph per fare la minimizzazione. Siccome dei nostri punti noi conosciamo anche l'errore, l'oggetto migliore per immagazzinare i dati e` un TGraphErrors, che e` completamente analogo a un TGraph, ma ha il metodo aggiuntivo
void SetPointError(int punto,double incertezzax,double incertezzay)
che permette di inserire degli errori su di un punto. Un frammento di codice per impacchettare i dati potrebbe essere:


  int Nlambda=7;
  double lambda[7]={579.1,577.0,546.1,491.6,435.8,407.7,404.7};
  double s=0.00055; // valore tipico per l'incertezza su n
  TGraphErrors g;
  for (int i=0; i<Nlambda; i++) {
    g.SetPoint(i,lambda[i],n[i]); // punto (lambda,n)
    g.SetPointError(i,0.,s);      // incertezza e' 0 su lambda e 0.00055 per n
  }

Passo 2: definizione della funzione che descrive i dati

Il passo successivo e' definire una funzione da adattare ai dati osservati. Questa funzione dovra` dare l'indice di rifrazione in funzione della lunghezza d'onda e dipendere dai parametri A e B.

ROOT usa la classe TF1 per descrivere una funzione parametrica monodimensionale.

Ci sono due possibilita` per definire una funzione TF1:

scriverne la formula
se la funzione e' semplice, possiamo semplicemente costruirla dandole un nome, l'espressione della funzione, e gli estremi di validita`. Nella formula si devono indicare con x la variable indipendente e con un numero tra parentesi quadre i vari parametri:
TF1 Cauchy("Cauchy","sqrt([0]+[1]/(x*x))",0.,1000.);
scrivere una funzione ausiliarla e promuoverla a TF1
se l'espressione della funzione e' piu` complicata, possiamo fare un passo intermedio e scrivere una funzione con la seguente interfaccia:
double funzione(double *x, double *p)
dove x indica il vettore delle variabili indipendenti (che per una funzione monodimensionale ha lunghezza 1) e p il vettore dei parametri (che nel nostro caso avra` lunghezza 2). Tale funzione puo` poi essere usata per costruire una TF1 dando un nome, un puntatore alla funzione, gli estremi di validita` ed infine il numero di parametri.
Nel nostro caso avremmo potuto farlo con il seguente frammento di codice:
```
double cauchyfun(double *x, double *p) {
  return sqrt(p[0]+p[1]/(x[0]*x[0]));
}
TF1 Cauchy("Cauchy",cauchyfun,0.,10000.,2)
```

Passo 3: condizioni iniziali

Al primo anno avete provato a fare fit usando delle funzioni predefinite (gaussiane e polinomi). Per queste funzioni speciali ROOT e' in grado di determinare un buon punto di partenza per la procedure di minimizzazione.

Per una funzione generica, questo non e` possibile e spetta all'utente fornire una stima iniziale dei parametri. Questo puo` venire fatto attraverso il metodo
void SetParameter(int parametro, double valore)
Siccome il valore di partenza alla fine puo` determinare il minimo che viene trovato, bisogna assicurarsi di darne una posizione verosimile. Una parte facoltativa dell'esercizio potrebbe essere verificare se il risultato della minimizzazione dipende o meno dai valori iniziali.

Nel nostro caso specifico, potremmo inizializzare i valori di A e B con il seguente frammento di codice:

Cauchy.SetParameter(0,2.);
Cauchy.SetParameter(1,50000.);

Passo 4: effettuare il fit

Una volta preparato tutto il necessario, effettuare la minimizzazione e` un'operazione banale: basta invocare il metodo:
int Fit(TF1* f1, char* option = "", char* goption = "", double xmin = 0, double xmax = 0)
dove i parametri hanno il seguente significato:

f1 puntatore alla funzione che descrive i dati nel grafico/istogramma;
option stringa contenente le opzioni del fit, di solito l'opzione piu` utile e la "Q" che dice di non stampare a schermo il risultato della minimizzazione. Se non si seleziona questa opzione, al termine del fit viene stampata una videata che indica per i vari parametri il valore ottenuto e l'incertezza stimata;
xmin e xmax definiscono l'intervallo in cui fare la minimizzazione, se questo non coincide con tutta l'estensione del grafico/istogramma.

Praticamente bastera` fare:
g.Fit(Cauchy,"Q")
semplicemente per evitare di avere un uscita intricata con stampe troppo prolisse.

Passo 5: estrarre i risultati

Una volta effettuato il fit, dobbiamo recuperare i risultati per stamparli su file nel formato che preferiamo o metterli in un istogramma.

Questo si puo` fare utilizzando i metodi
double GetParameter(int parametro)
double GetParError(int parametro)
della classe TF1:


A = Cauchy.GetParameter(0); // parametro A
sA= Cauchy.GetParError(0);  // errore sul parametro A
B = Cauchy.GetParameter(1); // parametro B
sB= Cauchy.GetParError(1);  // errore sul parametro B

Completare la simulazione

Modificare il programma che calcola gli indici di rifrazione, per completare ogni singolo esperimento con la determinazione di A e B. Questo si potrebbe realizzare costruendo una funzione separata:
void parametriCauchy(double& A, double& sA, double& B, double& sB, double* n)
che prende come input il vettore n degli indici di rifrazione misurati e restituisce i valori ottenuti nelle variabili A, sA, B e sB.

Tale funzione potete costruirla facilmente mettendo assieme nel modo opportuno tutti i frammenti di codice delle parti precedenti.

Anche in questo caso si tratta di valutare l'incertezza del metodo di misura dalla distribuzione degli scarti dei parametri A e B rispetto al proprio valore vero.

Si confronti il valore cosi' ottenuto con quello degli errori sA ed sB determinati da ROOT nei due casi:

tutte le sorgenti di errore sono attive,
solo l'errore su θ_m e` considerato.

Dovreste vedere che c'e' un palese disaccordo nel primo caso, mentre l'accordo e` quasi perfetto nel secondo.

N.B.: perche' i valori di sA ed sB e' necessario che il grafico venga costruito utilizzando le incertezze corretti su n(λ) e queste sono diverse per i due casi (nel secondo caso ci sono meno sorgenti di errore e quindi σ_n sara` piu` piccolo).

La motivazione di questo comportamento e' che, nel calcolo degli errori, ROOT assume che non ci sia correlazione tra i punti misurati. Questa assunzione e' vera nel secondo caso, ma non nel primo.

Il motivo per cui si assume assenza di correlazione tra i vari punti e' che le formule relative sono piu` semplici (anche le formule che il Taylor da' per i fit lineari sono valide solo in questa ipotesi).

E' possibile dare delle formule analitiche anche per il caso di punti correlati, ma sono cosi' scomode da usare che, alla fine, quasi sempre, una simulazione e' una strada piu` conveniente per ottenere la risposta corretta.

Esercizi

Per prendere dimestichezza con il fit una funzione parametrica a dei dati, provate a svolgere i seguenti esercizi.

Determinazione di parametri di una risonanza

Se prendiamo un oscillatore armonico forzato e smorzato, con equazione del moto:

l'ampiezza dell'oscillazione risultante e' data da:

In laboratorio sono state prese delle misure di ampiezza di oscillazione in funzione della frequenza di una forzante, con un'incertezza sulla singola misura di 0.01 unita` di ampiezza. Le trovate nel file datiRisonanza.txt, in cui la prima colonna indica le frequenza in Hz e la seconda le ampiezze in unita` arbitrarie.

Riordinare le righe del file in modo da mettere tutti i dati in ordine crescente di frequenza.
Fare un fit ai dati per determinare i parametri F, ω₀ e α del sistema (ricordate che la frequenza f vale ω/2π).
Confrontare i risultati ottenuti con i valori teorici dei paramettri:
- F = 2
- ω₀=11.31 rad/s
- α=0.03 s^-1

macro di soluzione

Simulatione della lettura di dati da un oscilloscopio

Supponiamo di avere un oscilloscopio che misura dei segnali di tensione con un'incertezza di 10 mV. Simulare una serie di misure dell'oscilloscopio se si osserva un segnale sinusoidale della forma:
V=V₀ cos( ωt+Φ₀)
dove:

V₀ = 1 V
la frequenza e` 0.5 Hz
Φ₀=0.6 rad
La misura consiste di 100 campionamenti distanziati di 0.1 s;

Presi i punti sperimentali cosi' ottenuti, provare a ricavare ampiezza, fase e frequenza facendo un fit ai dati: i valori ottenuti dal fit sono compatibili con i valori veri ottenuti dalla simulazione? Ripetere la procedura di misura e fit 100 volte e costruire un istogramma delle frequenze ottenute.

macro di soluzione

Relazione

Compilare la relazione