Introduzione Python MINISCUOLA DI FISICA DELLE PARTICELLE 2022

In questo breve esempio, senza la pretesa di imparare un linguaggio di programmazione in qualche ora, vediamo come possiamo in modo pratico utilizzare Python3 per maniploare, analizzare e visualizzare dati raccolti da rivelatori di particelle. Questi strumenti infatti non ci danno direttamente la quantità fisica di nostro interesse, tuttavia a partire da questi dati è possibile calcolarla o ricavarla

Download

• Questo quaderno è consultabile a questa pagina web
• Il quaderno eseguibile può essere scaricato da qui
• I dati utilizzati possono essere scaricati da qui

I quaderni jupyter

I quaderni jupyter sono un insieme di celle di markdown, dove si può scrivere del testo, eventualmente formattato (per esempio in corsivo o in grassetto), oppure delle formule matematiche $e^{i\pi}+1=0$, $\int_a^bf\left(x\right)\,dx = F\left(b\right)-F\left(a\right)$ e da celle di codice, dove le espressioni scritte vengono valutate

3+2

5

4*3

12

L'output, il risultato viene mostrato sotto. Ma solo dell'ultima espressione. Per stampare tutto quello che vogliamo, si può usare l'istruzione print, che stampa quello che si trova tra le sue parentesi tonde

4+5
3/2

1.5

print(4+5)
print(3/2)

9
1.5

Brevissima introduzione di Python

Le variabili

Le variabili sono dei contenitori che permettono di memorizzare un valore, sia esso un numero, una stringa (= una frase, del testo) o oggetti più complicati. Le variabili sono dunque costituite da un nome e un valore. Per creare una variabile, è sufficiente dare un nome e specificare il suo valore tramite l'operatore di assegnazione =.
Per le stringhe il valore deve essere racchiuso tra apici, singoli o doppi, mentre i numeri vengono scritti semplicemente

miaVariabile = 3
miaStringa = "Ciao"

Per vedere cosa c'è dentro una stringa

print(miaVariabile)
print(miaStringa)

3
Ciao

Per modificare il valore è sufficiente effettuare una nuova assegnazione

miaStringa = "Buongiorno"
print(miaStringa)

Buongiorno

Le variabili sono utili per contenere il risultato di certe operazioni. Si possono fare le usuali operazioni aritmetiche: +, -, *, /. Per l'elevamento a potenza, si deve utilizzare un doppio asterisco **.

Sono validi tutti i nomi alfanumerici, senza spazi, e i trattini bassi (underscore). Il nome non può però iniziare con un numero. Python è case-sensitive, ovvero le lettere maiuscole sono importanti

È buona norma dare alle variabili dei nomi esplicativi, che permettano di capire cosa c'è dentro (evitando per esempio a,b come fatto in questo esempio didattico)

Ci sono due modalità principali: separare con un underscore le parole oppure utilizzare la notazione a cammello

a = 3
b = 7

somma = a+b
prodotto = a*b
potenza = 3**7

print(somma)
print(prodotto)
print(potenza)

Somma = 29

# Con la virgola posso stampare due output
# Il cancelletto in una riga di codice identifica un COMMENTO, ovvero
# una nostra annotazione, non viene eseguito
print(Somma, somma)

10
21
2187
29 10

Le f-strings

Si tratta di un modo comodo per stampare una stringa formattata. Anteponendo una f davanti agli apici, tutto quello che si mette dentro parentesi graffe viene valutato, non stampato letteralmente. Riprendiamo l'esempio di prima

a = 3
b = 7

somma = a+b
prodotto = a*b
potenza = 3**7

# Si noti la f davanti agli apici
print(f"La somma di {a} con {b} vale {somma}, il prodotto vale {prodotto} e la potenza vale {potenza}")

La somma di 3 con 7 vale 10, il prodotto vale 21 e la potenza vale 2187

Nelle stringhe sono ammessi dei caratteri speciali, detti caratteri di escape: \t rappresenta una tabulazione, \n per andare a capo. Nel caso in cui non si volesse tale comportamento, si antepone una r agli apici.

Le f-strigs sono molto potenti. È possibile per esempio controllare il numero di decimali che vengono visualizzati

print(f"7/3:\t{7/3}\t{7/3:.4f}\n5/9:\t{5/9}\t{5/9:.2f}")

7/3:	2.3333333333333335	2.3333
5/9:	0.5555555555555556	0.56

Le liste

Le liste sono da pensare come una collezione di variabili, tutte dello stesso tipo, ma che anziché avere un nome ciascuna, si dà solamente un nome all'insieme. Gli elementi sono ordinati, e infatti si accede ad un valore mediante l'indice. Le liste sono racchiuse da parentesi quadre. L'indice in qualsiasi contesto è sempre racchiuso tra parentesi quadre.

ATTENZIONE: In python, come nella maggior parte dei linguaggi di programmazione, la prima posizione corrisponde all'indice zero, la seconda all'indice 1 e via dicendo

Esiste la possibilità di utilizzare indici negativi per iniziare a contare dal fondo. In questo caso -1 rappresenta l'ultimo elemento

Colgo l'occasione per sottolineare come ho utilizzato un apostrofo nella f-string, che sarebbe l'apice singolo. La flessibilità nel poter usare sia gli apici singoli che i doppi apici nasce proprio da questo: essendo nel mio caso i doppi apici il delimitatore, l'apice singolo è un nome valido. Per esigenze più particolari, con l'operatore + si possono concatenare due stringhe tra di loro

# Creo una lista di numeri
listaNum = [3,4,7,5,1]

# Stampo la lista
print(listaNum)

# Stampo il secondo elemento (indice 1)
# Accedo all'elemento con le parentesi quadre
print(listaNum[1])


# Stampo la somma del primo e dell'ultimo elemento
print(f"Somma del primo e dell'ultimo elemento: {listaNum[0] + listaNum[-1]}")

[3, 4, 7, 5, 1]
4
Somma del primo e dell'ultimo elemento: 4

Ancora sulle liste:

• Con il metodo len posso ottenere la lunghezza di una lista
• Con il metodo append posso appendere degli elementi alla fine

Si possono fare molte altre cose, tipo inserire in una posizione specifica, rimuovere in una posizione specifica, ma queste sono le operazioni più usuali

# Posso ottenere il numero di elementi di una lista
print(len(listaNum))

# Aggiungo un numero alla lista
listaNum.append(5)
print(listaNum)

listaNum.append(27)
print(listaNum)

# Elimino ogni volta l'ultimo valore
listaNum.pop()
print(listaNum)
listaNum.pop()
print(listaNum)
listaNum.pop()
print(listaNum)

5
[3, 4, 7, 5, 1, 5]
[3, 4, 7, 5, 1, 5, 27]
[3, 4, 7, 5, 1, 5]
[3, 4, 7, 5, 1]
[3, 4, 7, 5]

Si può anche selezionare un intervallo mediante inizio:fine, con la convenzione che il punto finale si intende escluso

listaNum.append(5)
listaNum.append(27)
print(listaNum)

# Selezioni gli indici 3 e 4
print(listaNum [3:5])

[3, 4, 7, 5, 5, 27]
[5, 5]

Posso anche fare una lista di stringhe

lstStr = ["ciao", "come", "stai?"]
print(lstStr)

# Uso il + per concatenare assieme tutte le stringhe

print(lstStr[0] + ", " + lstStr[1] + " " + lstStr[2])

['ciao', 'come', 'stai?']
ciao, come stai?

I numeri decimali

Finora abbiamo visto le stringhe e i numeri interi. Python è in grado di gestire i numeri decimali (floating point) ed anche la notazione scientifica (la e sta per il numero di zero da mettere dopo)

a = 3.0
b = 7.5
print(f"{a+b}")

# Anche mischiando tipi di numeri, python se la cava egregiamente
a = 7
b = 12.34
print(f"{a+b}")


velocitaLuce = 2.97e8 #m/s
print(velocitaLuce)

10.5
19.34
297000000.0

Le strutture decisionali

Finora abbiamo eseguito semplici assegnazioni e/o operazioni aritmetiche. Noi vogliamo che il nostro programma sappia prendere delle decisioni. Le strutture if ... else ci permettono di creare una biforcazione per il nostro programma, e verrà scelta solo una delle due strade, a seconda che una data condizione sia vera o falsa.

È il momento buono per introdurre un altro tipo di dati, i boolean: si tratta di dati che possono essere solamente True oppure False (notare l'iniziale maiuscola). Su questi dati esistono alcune operazioni binarie (in quanto riguardano il confronto di due variabili) standard come AND (&) (vero solo se entrambi sono veri) oppure OR (|) (vero solo se almeno uno dei due è vero). L'operazione unaria è la negazione NOT

bool1 = True
bool2 = False

print(f"AND: {bool1 & bool2}\nOR:{bool1 | bool2}")
print(f"NOT: {not bool1}\nOR:{not bool2}")

AND: False
OR:True
NOT: False
OR:True

È possibile ottenere un risultato booleano confrontando due variabili, in particolare testando la loro uguaglianza == (singolo uguale per assegnare un valore, doppio uguale per testare l'uguaglianza), oppure >, >=, <, <=

a = 3
b = 4

print(a>3)
print(a==a)
print(a==b)

False
True
False

if ... else

Vediamo quindi come far prendere al nostro programma una diversa strada a seconda che una certa condizione sia vera oppure falsa
Si noti l'importanza dell'indentazione, ovvero il fatto che il blocco di codice sotto l'if sia spostato. Questa cosa è fondamentale, pena un errore / il non corretto funzionamento del nostro programma

a = 3
#a = 13
b = 7

if a>b:
    print("è maggiore")
else:
    print("è minore")

è minore

elif (else if)

Tramite questa keyword è possibile verificare più condizioni in sequenza

a = 3
#a = 13
b = a

if a>b:
    print("è maggiore")
elif a<b:
    print("è minore")
else:
    print("è uguale")

è uguale

I cicli for

Supponiamo di avre una lista, che contiene tanti elementi, magari i dati raccolti da un rivelatore (sviariate migliaia) di particelle e vogliamo chiamarli in causa tutti, ma senza dover scrivere una riga per ciascun elemento, vogliamo fare in modo che si chiamino automaticamente...

I cicli for ci permettono di eseguire una certa operazione per tutti gli elementi di un oggetto iterabile, i seguenti esempi spero chiariranno l'idea

listaLunga = ["a", "b", "c", "d", "e", "f", "g", "h", "i"] # ...

for lettera in listaLunga:
    print(lettera)

a
b
c
d
e
f
g
h
i

Per avere un contatore di numeri si può usare la funzione range(x). Si noti che genera $x$ numeri, da 0 a $x-1$

for i in range(10):
    print(f"Il quadrato di {i} vale {i**2}")

Il quadrato di 0 vale 0
Il quadrato di 1 vale 1
Il quadrato di 2 vale 4
Il quadrato di 3 vale 9
Il quadrato di 4 vale 16
Il quadrato di 5 vale 25
Il quadrato di 6 vale 36
Il quadrato di 7 vale 49
Il quadrato di 8 vale 64
Il quadrato di 9 vale 81

I cicli while

Un'altra tipologia di ciclo sono i cicli while. A differenza dei cicli for, dove si va ad iterare sugli elementi di un certo insieme, il ciclo while continua a ripetersi finche una certa condizione è vera, ma non si sa dire a priori quante volte si dovrà ripetere.

Attenzione a fare in modo che la condizione prima o poi diventi falsa, altrimenti si avrà un ciclo infinito da cui non si esce mai, e rischia di far bloccare il PC

numero = 54382832 # Un numero grande

# Lo dimezzo finchè non ottendo un numero minore di 1
# non so dire a priori quante volte si deve ripetere
while numero>1:
    numero  = numero /2

print(numero)

0.8103673458099365

Alcuni esempi avanzati, ma che possono tornare utili

List comprehension

Supponiamo di voler creare una lista che contenga i quadrati dei numeri da 0 a 9. A partire dall'esempio di prima, uno potrebbe fare

# Inizializzo una lista vuota
lstQuadrati = []

for i in range(10):
    lstQuadrati.append(i**2)
    
print(lstQuadrati)

[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

Le List comprehension sono una notazione compatta per fare lo stesso lavoro (in realtà sono anche più efficienti...)

lstQuadratiEff = [i**2 for i in range(10)]
print(lstQuadratiEff)

[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]

Zip

Supponiamo di avere due liste e di voler iterare sulle coppie, avendo ad ogni iterazione la coppia in posizione 0, poi la coppia in posizione 1... Si potrebbe pensare di ciclare sugli indici

lst1 = ["a", "b", "c", "d", "e", "f", "g", "h", "i"]
lst2 = ["A", "B", "C", "D", "E", "F", "G", "H", "I"]

for i in range(len(lst1)): # Questa è una delle ragioni per cui l'ultimo numero è escluso
    print(f"{lst1[i]}\t{lst2[i]}")

a	A
b	B
c	C
d	D
e	E
f	F
g	G
h	H
i	I

Utilizzando il metodo zip si può ottenere in modo leggermente più comodo lo stesso effetto. Eliminando un indice, si elimina un elemento di complessità, e il codice nel suo inisieme risulterà più leggibile

for lett1, lett2 in zip(lst1, lst2):
    print(f"{lett1}\t{lett2}")

a	A
b	B
c	C
d	D
e	E
f	F
g	G
h	H
i	I

Enumerate

Supponiamo infine di voler ciclare su una lista, ma di voler anche sapere in che indice siamo. Nuovamente, si potrebbe pensare di ciclare sugli indici della lista nel seguente modo

lst1 = ["a", "b", "c", "d", "e", "f", "g", "h", "i"]

for i in range(len(lst1)):
    print(f"In posizione{i} si trova l'elemento {lst1[i]}")

In posizione0 si trova l'elemento a
In posizione1 si trova l'elemento b
In posizione2 si trova l'elemento c
In posizione3 si trova l'elemento d
In posizione4 si trova l'elemento e
In posizione5 si trova l'elemento f
In posizione6 si trova l'elemento g
In posizione7 si trova l'elemento h
In posizione8 si trova l'elemento i

Il metodo enumerate ci permette di ottenere lo stesso risultato in un modo leggermente più pratico

lst1 = ["a", "b", "c", "d", "e", "f", "g", "h", "i"]

for indice, valore in enumerate(lst1):
    print(f"In posizione{indice} si trova l'elemento {valore}")

In posizione0 si trova l'elemento a
In posizione1 si trova l'elemento b
In posizione2 si trova l'elemento c
In posizione3 si trova l'elemento d
In posizione4 si trova l'elemento e
In posizione5 si trova l'elemento f
In posizione6 si trova l'elemento g
In posizione7 si trova l'elemento h
In posizione8 si trova l'elemento i

Piccolo break

A questo punto abbiamo visto i concetti fondamentali di Python. Sebbene in linea di principio con questi strumenti possiamo fare tutto, esistono dei moduli ulteriori che permettono di estendere le funzionalità. Di seguito vedremo brevemente come utilizzare i tre metodi scientifici principali:

• numpy: per gestire in modo efficiente liste di numeri
• matplotlib: per creare grafici, per visualizzare i dati
• scipy: per alcuni strumenti matematici avanzati (fit) che ci torneranno utili

Breve introduzione di numpy

Prima di tutto i moduli esterni devono essere importati mediante la keyword import. Nel caso di numpy è sufficiente scriverer import numpy. Si può anche dare un alias (una sorta di soprannome) alla libreria importata. L'alias lo si può scegliere a proprio piacimento, ma è convenzione universalmente accettata che numpy venga importato come np.

# Importo il modulo
import numpy as np

L'oggetto principale di numpy sono una sorta di liste più potenti, detti ndarray. Si creano a partire da una lista, chiamando np.array()

arr = np.array([1,2,3,4,5,6,7])
print(arr) # Mi stampa i valori
arr # Output di default, mi dice anche che è un array

[1 2 3 4 5 6 7]

array([1, 2, 3, 4, 5, 6, 7])

Gli array hanno alcuni campi che mi possono dare alcune informazioni

print(arr.shape, arr.dtype)

(7,) int32

Il loro punto di forza risiede nellle operazioni vettorializzate: è possibile effettuare operazioni elemento ad elemento, senza dover fare dei cicli for. Oltre alla comodità di scrittura, ha delle forti conseguenze in termini di performance (sono molto più veloci). Vediamo alcuni esempi

# Sommo un certo numero a tutti gli elementi
print(arr+3)

# Elevo al quadrato
print(arr**2)

# Ne faccio la radice quadrata
print(arr**.5)

[ 4  5  6  7  8  9 10]
[ 1  4  9 16 25 36 49]
[1.         1.41421356 1.73205081 2.         2.23606798 2.44948974
 2.64575131]

Alcune funzioni

Esistono alcuni metodi (alcune funzioni) privilegiati che possono essere chiamate direttamente dagli array, vediamo alcuni esempi

print(arr.max()) # Massimo dell'array
print(arr.min()) # Minimo dell'array

print(arr.mean()) # Valor medio
print(arr.std()) # Deviazione standard

7
1
4.0
2.0

Tutte queste funzioni, più molte altre, possono essere chiamate con la notazione usuale

# Funzioni che ritornano un numero a partire da un vettore
print(np.max(arr))
print(np.min(arr))
print(np.mean(arr))
print(np.std(arr))

print(np.median(arr)) # Mediana


# Funzioni che agiscono elemento per elemento, ritornando un vettore
print(np.sqrt(arr)) # Radice quadrata (elemento per elemento)

# Esistono per esempio tutte le funzioni trigonometiche: sin, cos, tan...

7
1
4.0
2.0
4.0
[1.         1.41421356 1.73205081 2.         2.23606798 2.44948974
 2.64575131]

Esistono anche altre funzioni comode per generare vettori di numeri

arange

# Vettore di 100 punti (ultimo escluso)
np.arange(100)

array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
       17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33,
       34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50,
       51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67,
       68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
       85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99])

# Vettore da 30 a 100 a step di 3
np.arange(30, 100, 3)

array([30, 33, 36, 39, 42, 45, 48, 51, 54, 57, 60, 63, 66, 69, 72, 75, 78,
       81, 84, 87, 90, 93, 96, 99])

linspace

# 100 punti equispaziati tra 0 e 10 
# In questo caso l'estremo superiore è incluso, ma si può 
# escludere con un certo parametro opzionale
np.linspace(0, 10, 100)

array([ 0.        ,  0.1010101 ,  0.2020202 ,  0.3030303 ,  0.4040404 ,
        0.50505051,  0.60606061,  0.70707071,  0.80808081,  0.90909091,
        1.01010101,  1.11111111,  1.21212121,  1.31313131,  1.41414141,
        1.51515152,  1.61616162,  1.71717172,  1.81818182,  1.91919192,
        2.02020202,  2.12121212,  2.22222222,  2.32323232,  2.42424242,
        2.52525253,  2.62626263,  2.72727273,  2.82828283,  2.92929293,
        3.03030303,  3.13131313,  3.23232323,  3.33333333,  3.43434343,
        3.53535354,  3.63636364,  3.73737374,  3.83838384,  3.93939394,
        4.04040404,  4.14141414,  4.24242424,  4.34343434,  4.44444444,
        4.54545455,  4.64646465,  4.74747475,  4.84848485,  4.94949495,
        5.05050505,  5.15151515,  5.25252525,  5.35353535,  5.45454545,
        5.55555556,  5.65656566,  5.75757576,  5.85858586,  5.95959596,
        6.06060606,  6.16161616,  6.26262626,  6.36363636,  6.46464646,
        6.56565657,  6.66666667,  6.76767677,  6.86868687,  6.96969697,
        7.07070707,  7.17171717,  7.27272727,  7.37373737,  7.47474747,
        7.57575758,  7.67676768,  7.77777778,  7.87878788,  7.97979798,
        8.08080808,  8.18181818,  8.28282828,  8.38383838,  8.48484848,
        8.58585859,  8.68686869,  8.78787879,  8.88888889,  8.98989899,
        9.09090909,  9.19191919,  9.29292929,  9.39393939,  9.49494949,
        9.5959596 ,  9.6969697 ,  9.7979798 ,  9.8989899 , 10.        ])

Caricare dei dati

Il modo in cui un fisico usa Python è per effettuare l'analisi di dati proveniente da qualche strumento scientifico, nel nostro caso saranno rivelatori di particelle. La maggior parte dei dati sono in formato ASCII, organizzati in colonne, dove ogni colonna contiene una certa quantità fisica, mentre ogni riga rappresenta un evento differente.

Il file run020252.dat ha la seguete forma

001618564101 91 1297 67 107 
001618564101 378 82 66 106 
001618564101 364 90 67 107 
001618564101 370 78 66 106 
001618564101 87 1034 66 107 
001618564101 92 222 67 107 
001618564101 90 1162 67 107 
001618564101 568 80 67 107 
001618564101 823 84 68 107 
001618564101 627 76 67 107 
001618564101 92 669 66 107 
001618564101 1323 80 67 106 
001618564101 2079 85 68 107 
001618564101 93 1627 67 108 
001618564101 90 537 67 107 
001618564101 1320 80 67 107 
001618564101 307 76 66 107 
001618564101 553 74 67 107 
...

Per caricarlo è sufficiente usare il metodo np.loadtxt(). Ci ritorna una matrice che contiene tutti i dati

dati = np.loadtxt("run020252.dat")
print(dati.dtype, dati.shape)

float64 (300001, 5)

Per ottenere la secondoa e la terza colonna, si potrebbe fare

a = dati[:,1] # Tutte le righe, colonna 1
b = dati[:,2] # Tutte le righe, colonna 2

matriceCol1e2 = dati[:,1:3] # Ricordiamo che con i :, l'ultimo indice è escluso

print(a.shape, b.shape, matriceCol1e2.shape)

(300001,) (300001,) (300001, 2)

Per caricare ogni colonna in un vettore diverso basta specificare unpack = True

_, pur essendo un nome legittimo per una variabile, viene convenzionalmente utilizzato per una variabile che non ha alcuna importanza. Infatti se voglio effettuare l'unpack dei dati, sono obbligato a fornire tante variabili quante sono le colonne del mio file

_, col1, col2, _, _ = np.loadtxt("run020252.dat", unpack = True)
print(col1.shape, col2.shape)

(300001,) (300001,)

Breve introduzione di matplotlib

Questo modulo, in particolare il sottomodulo pyplot consente di creare grafici a partire da vettori. Iniziamo importando il modulo. Anche qui esiste un alias universalmente condiviso

import matplotlib.pyplot as plt

Alcuni esempi di come creare un grafico

Iniziamo a crare dei vettori, degli insiemi di $x$ e $y$ da rappresentare

xVect = np.linspace(0, 2*np.pi, 1000) 

ySin = np.sin(xVect)
yCos = np.cos(xVect)

yStrano = np.sin(xVect + np.pi/4 )

Esempio super-minimo (ma da evitare)

plt.plot(xVect, ySin)

[<matplotlib.lines.Line2D at 0x1d155015520>]

Esempio minimo

fig, ax = plt.subplots()
ax.plot(xVect, ySin)
plt.show()

Vediamo cos'altro si può fare

Cerco di fare un esempio completo commentando nel codice le varie personalizzazioni possibili
Questo si tratta di un esempio. Consultando le documentazioni online si vede come per ciascuna funzione esistano molti più parametri impostabili

Tranne la prima e l'ultima riga, tutti i comandi possono essere messi nell'ordine che si vuole (Unico caveat, quando si chiama il comando legend viene stampato solo quello che è stato chiamato sopra)

# Creo la figura in cui andrò a disegnare (prendetelo per vero)
fig, ax = plt.subplots()

# Dimensioni (in pollici)
fig.set_size_inches(12,5)

# Disegno una curva
# Specifico il vettore delle x, il vettore delle y, il tipo di linea, il colore, la legenda
ax.plot(xVect, ySin, linestyle = "-.", color = "tab:red", linewidth = 3, label = "$y=\sin(x)$")

# Per disegnare un'altra curva, basta richiamare il metodo plot
# Per i parametri sopra, esistono anche degli alias compatti
ax.plot(xVect, yCos, ls = "--", c = "tab:green", lw = 2, label = "$y=\cos(x)$")

# Inoltre, per il marker, il colore e il tipo di linea, esiste una notazione compatta
ax.plot(xVect, yStrano, ":k")



# Etichette assi
ax.set_xlabel("Ascissa", fontsize = 14)
ax.set_ylabel("Ordinata", fontsize = 14)

# Titolo
ax.set_title("Plot di funzioni trigonometriche", fontsize = 16, color = "red", fontweight = "bold")


# Griglia
ax.grid()

# Creo la legenda, specificando titolo e dimensioni
ax.legend(title = "Legenda", fontsize = 14)



# Mostro la figura con tutte le cose che ci ho disegnato
plt.show()

Un esempio di due plot in una singola figura

fig, ax = plt.subplots(2,1) # 2 righe, una colonna

fig.set_size_inches(12,5)

# Permette di gestire la spaziatura tra due grafici di una stessa figura
fig.subplots_adjust(hspace = .4)

#questa volta abbiamo una sola figura, ma ax è un vettore di assi, uno per ciascun subplot

ax[0].plot(xVect, ySin, c = "tab:green", label = "$y = \sin(x)$")
ax[1].plot(xVect, yCos, c = "tab:green", label = "$y = \cos(x)$")

# Essendo un vettore ax, per impostare cose comuni a ciascun asse, posso iterare

for a in ax:
    a.grid()
    
    a.set_xlabel("Ascissa", fontsize = 14)
    a.set_ylabel("Ordinata", fontsize = 14)
    
    a.legend(title = "Legenda", fontsize = 14)

plt.show()

Appendice: Mettiamoci all'opera

Di seguito si mostra rapidamente come vengono gestiti in python due degli oggetti fondamentali con cui un fisico avrà a che fare: istogrammi e fit.

Gli istogrammi permettono di studiare la distrubizione di un set di dati: vengono fissati dei bordi, ogni coppia contigua di bordi definisce un bin. L'istogramma rappresenta, per ciascun bin, il numero di eventi che si trovavano in quel dato intervallo

Il fit dei dati consiste nel fittare un certo set di dati oppure un certo istogramma con una data funzione matematica (retta, parabola, esponenziale, gaussiana...) i cui parametri hanno generalmente un significato fisico.

Istogramma e fit gaussiano

In questo esempio didattico, che commenterò durante il codice,

1. Inizio generando un vettore di punti distribuiti in modo gaussiano (la distribuzione gaussiana è conosciuta anche come distribuzione normale)
2. Procedo alla costruzione dell'istogramma
3. Effettuo il fit dell'istogramma con una funzione gaussiana definita da me

Step 1: generazione numeri

Il modulo np.random contiene una serie di funzioni utili per la generazione di numeri casuali. La distribuzione normale, in particolare, dipende da due parametri: il valor medio $\mu$ e la deviazione standard $\sigma$. (Hint: scrivere np.random. e premere tab per usare la funzione di auto-completamendo, Shift + tab per visualizzare la documentazione)

Come potete vedere, generare un milione di punti è facilissimo

# Genero un vettore di 1000000 numeri casuali con mu = 3 e sigma = 5
vectGauss = np.random.normal(3, 5, 1000000)

Step 2: istogramma

Per costruire un istogramma si può utilizzare la funzione np.histogram. All'ordine zero, chiede che gli si passi il set di dati da istogrammare. Ulteriori parametri configurabili sono il numero di bin o l'intervallo entro cui considerare i punti (quelli fuori vengono scartati). Questa funzione ritorna il numero di conteggi per ciascun bin e i bordi (per ottenere il centro serve uno step ulteriore)

Successivamente vado a plottare l'istogramma da noi costruito: vediamo che ha proprio la forma di una gaussiana centrata sul 3

# Costruisco l'istogramma specificando il numero di bin. 
# Volendo, uno potrebbe specificare anche l'intervallo da considerare
h, bins = np.histogram(vectGauss, bins = 100)#, range = (-7, 13))

# Trovo il centro dei bin, spostando tutti i bordi sinistri di metà ampiezza
binc = bins[:-1] + (bins[1]-bins[0]) / 2






# Plotto l'istogramma
fig, ax = plt.subplots()
fig.set_size_inches(12,5)


# Si noti il parametro ds = "steps-mid" (drawstyle) che specifica di non congiungere
# i punti con una linea retta, ma di disegnare l'istogramma a barre
ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", label = "Istogramma" )
        
        
# Si tratta di una funzione carina che permette di colorare la regione sottesa ad una certa curva
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)


ax.set_title("Spettro", fontsize = 16)
ax.set_xlabel("Energia [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)

ax.grid(True)
 
# Questo parametro ora non ci serve, ma è il modo in cui si imposta l'asse y logaritmico
# Ne parleremo quando ci servirà
#ax.set_yscale("log")

# Un altro paio di opzioni per la legenda, ad esempio personalizzando la dimensione del titolo
# e scegliendo dove posizionarla
ax.legend(fontsize = 14, title = "Legenda", loc = "upper right", title_fontproperties = {"size":16})

plt.show()

Step 3: fit dei dati

Come anticipato sopra, supponiamo che quei dati siano stati raccolti da un rivelatore di particelle, e che il valor medio abbia un significato fisico. Noi vogliamo trovare quale sia la gaussiana (nel senso di quali siano i parametri) che meglio approssimi la nostra distribuzione.

Iniziamo ricordando che la funzione Gaussiana ha la seguente espressione matematica \begin{equation} G(x) = A * e^{\dfrac{-\left(x-\mu\right)^2}{2\sigma^2}} \end{equation} dove $A$ è una costante di normalizzazione, $mu$ il valor medio e $\sigma$ la deviazione standard.

Senza la pretesa di fare una trattazione rigorosa, descrivo come operativamente si possa effettuare il fit

• Si definisce una funzione di python (così imparate come si fa) che, dati come argomenti la variabile $x$ (al primo posto) e successivamente tutti i parametri, ritorni il valore della funzione. Per definire una funzione serve la keyword def, quindi tra parentesi si mettono i parametri ed indentato si scrive il codice che la definisce. La funzione fa una serie di operazioni basate sui parametri di input e ritorna il valore definito dalla keyword return

• Per effettuare i fit ci servirà il modulo scipy, in particolare la funzione scipy.optimize.curve_fit

• Effettuo il fit chiamando la funzione e specificando i dati e la funzione da fittare. Mi verrà restituito un vettore con i parametri della gaussiana che meglio approssimano la mia funzione. Un yema molto importante in statistica sono gli errori, nel senso che bisogna specificare ciascun punto che vado a fittare quanto pesa, quanto è importante. Ci sono varie considerazioni statistiche, ma per i dati provenienti da rivelatori di particelle, in generale (o almeno in prima approssimazione) si attribuisce come incertezza sulla misura la radice del numero stesso

• Vado a disegnare lo stesso plot di prima a cui sovrappondo anche la curva fittata

# Con i tre apici si fanno stringhe su più righe
# è un modo pratico per fare commenti lunghi nel codice
"""
Definisco una funzione che prende in input la variabile x e poi 
tutti i parametri. Essa ritornerà il valore della funzione
"""
def myGauss(x, a, mu, sigma):
    return a * np.exp( (-(x-mu)**2) / (2*sigma**2) )


# Importo una data funzione da un certo modulo
from scipy.optimize import curve_fit


"""
Effettuo il fit specificando nell'ordine
- La funzione da fittare
- I dati x e y
- Le incerteze sulle y
- Un parametro da "prendere per vero", serve a dire che vogliamo 
davvero usare quelle incertezze
- p0 serve a fornire una stima ragionevole sul punto di partenza
Per distribuzioni complicate, senza un aiutino, è complicato far
convergere il fit.
- maxfev serve a dire quante iterazioni fare

In questo caso viene inizializzato in modo smart... argmax ritorna 
l'indice in corrispondenza di cui c'è il massimo


Oggetti ritornati
- popt è il vettore dei parametri ottimizzati
- pcov è la "matrice delle covarianze", un oggetto complicato, ma i 
cui elementi diagonali sono i quadrati degli errori associati ai 
parametri estratti dal fit
"""
p0 = [np.max(h), binc[np.argmax(h)], 5]
# popt, pcov = curve_fit(myGauss, binc, h, sigma = np.sqrt(h), 
#                       absolute_sigma = True, 
#                       p0 = p0, maxfev = 10000000)


"""
Lascio l'esempio perchè è giusto, ma in realtà in questo modo si darebbe
molta più importanza alle code e il fit non convergerebbe. In questo caso
quel tipo di errore non è giustificato...
"""
popt, pcov = curve_fit(myGauss, binc, h, p0 = p0)


# Stampo i parametri con i relativi errori
print(f"Il parametro A è risultato essere {popt[0]} +- {np.sqrt(pcov[0,0])}")
print(f"Il parametro mu è risultato essere {popt[1]} +- {np.sqrt(pcov[1,1])}")
print(f"Il parametro sigma è risultato essere {popt[2]} +- {np.sqrt(pcov[2,2])}")



# Plot dei dati di prima a cui sovrappongo la gaussiana fittata
fig, ax = plt.subplots()
fig.set_size_inches(12,5)


ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", label = "Istogramma" )
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)



"""
In questa riga vado a disegnare la gaussiana fittata. In particolare, come x uso le stesse x, mentre
le y le calcolo a partire dal mio vettore di x e i parametri estratti dal fit.
L'asterisco davanti ad un vettore serve per effettuare l'unpacking, è come se scrivessi 
popt[0], popt[1], popt[2] (e viene fatto in automatico per tutti gli elementi del vettore)
Infatti la funzione myGauss vuole 4 argomenti: le x e i tre parametri A, mu, sigma
"""
ax.plot(binc, myGauss(binc, *p0), ls = "--", c = "k", lw = 3, label = "Gauss fittata")


ax.set_title("Spettro", fontsize = 16)
ax.set_xlabel("Energia [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)

ax.grid(True)
ax.legend(fontsize = 14, title = "Legenda", loc = "upper right", title_fontproperties = {"size":16})

plt.show()

print (p0)

Il parametro A è risultato essere 37145.594169789794 +- 23.606783705955472
Il parametro mu è risultato essere 2.995679403330614 +- 0.0036677073032359746
Il parametro sigma è risultato essere 4.997992364966411 +- 0.0036677074222852836

[37320, 2.8607893596105196, 5]

Addendum

È bene notare come la parte colorata di verde chiaro sia solo estetica, i dati sono rappresentati dalle altezze di ciascun bin. Per chiarire meglio rifaccio lo stesso plot come uno scatter plot, ovvero disegnando dei singoli punti senza linee di congiungimento. Sono questi punti che sono stati fittati con la gaussiana.

fig, ax = plt.subplots()
fig.set_size_inches(12,5)


# Con ls = "" specifico che non voglio la linea, mentre marker specifica
# il tipo di marker. (ms = marker size...)
ax.plot(binc, h,  c = "tab:green", label = "Istogramma",
       ls = "", marker = "*")#, ms = 10)


ax.set_title("Spettro", fontsize = 16)
ax.set_xlabel("Energia [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)

ax.grid(True)
ax.legend(fontsize = 14, title = "Legenda", loc = "upper right", title_fontproperties = {"size":16})

plt.show()

Conclusioni

Cosa ho omesso

• (Funzioni in Python)
• Blocchi swith case
• Leggere/salvare npy/npz