Lezione 5

Introduzione

Scopo di questa sessione di laboratorio è di 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.

Il processo logico che si segue nell'impostare una simulazione è 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 quantità, 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 così calcolati sono una valutazione corretta dei valori veri delle costanti fisiche?
qual è 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 può 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.

**Figura 1:** goniometro usato nell'esperimento dello spettrometro a prisma.

L'apparato è indicato in figura 1 e consiste in un goniometro sul quale viene posizionato il prisma. Una lampada a vapori di mercurio viene posizioneta da un lato del canocchiale con due collimatori per produrre un fascio luminoso che incide sul prisma. Il fascio di luce riflesso o rifratto viene osservato tramite un altro canocchiale. Gli angoli corrispondeni all'orientamento del supporto sul goniometro ed alla posizione dei canocchiali sono leggibili su di un nonio.

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 già effettuato l'esperimento in laboratorio, può 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, secondo la procedura illustrata in figura 2.
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.
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:
α= 180^o-(θ₁-θ₂)

Figura 2: Misura dell'angolo di apertura del prisma attraverso la riflessione.

$\resizebox*{1\columnwidth}{!}{\includegraphics{prisma.eps}}$
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 λ [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

**Figura 2:** 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.

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 così 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.: perché i valori di sA ed sB siano stimati correttamente, è necessario che il grafico venga costruito utilizzando le incertezze corrette su n(λ) e queste sono diverse per i due casi (nel secondo caso ci sono meno sorgenti di errore e quindi σ_n sarà più piccola).

La motivazione di questo comportamento è che, nel calcolo degli errori, ROOT assume che non ci sia correlazione tra i punti misurati. Questa assunzione è 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 così scomode da usare che, alla fine, quasi sempre, una simulazione e' una strada piu` conveniente per ottenere la risposta corretta.

Esercizi sui fit ai dati

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 veri dei parametri utilizzati per creare l'insieme delle misure:
- F = 2
- ω₀=11.31 rad/s
- α=0.03 s^-1

macro di soluzione

Per avere un'idea se parametri stimati sono correlati o meno, si può provare a fissarne uno alla volta, usando il metodo FitParameter, invece di SetParameter usato in precedenza. In caso di correlazione, il fissare uno dei parametri al valore vero, risulta in uno spostamento anche degli altri.

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;
Ad ogni misura bisogna assegnare un errore estratto da una \ distribuzione gaussiana con σ=0.01 V.

Presi i punti sperimentali così 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, da confrontare poi con il valore dell'incertezza stimato da RooT.