Assorbimento Pb

marco "operativa"

Apre tutti i file, legge gli adc necessari da tutti i file .dat e .log2 (necessario per avere il tempo di acquisizione), calcola il numero di eventi sotto il picco (fit gaussiano)

In [1]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt




### write a function to evaluate deadtime
def compute_corr_factor(nrun):
    datafile = './dati/run%6.6d.log2' % nrun
    time1, time2, time3, time4 = np.loadtxt(datafile, unpack=True)
    # time3 contiene i conteggi veri, dividi per i "presunti" e ottieni il fattore di correzione
    corr = sum(time3)/20000
    return corr

## gaussiana
def gaus(x, A, mu, sigma):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

### write a function to read data and compute the event rate
def read_data_and_compute_rate(nrun):
    datafile = './dati/run%6.6d.dat' % nrun
    time, adc1, adc2, adc3, i5 = np.loadtxt(datafile, unpack=True)
    elapsed_time = time[-1]-time[0]
    hrange=(1000, 2500) # <-- it's a tuple
    fig, ax = plt.subplots(1)
    plt.title('Run %d'%nrun, fontsize=18)
    plt.xlabel('ADC', fontsize=16)
    plt.ylabel('number of events', fontsize=16)
    h, bins, pathes = plt.hist(adc1, 500, histtype='step', range=hrange, edgecolor='black', linewidth=1.2, label='data')
    # find subrange for fitting
    subrange=(1670, 1850)
    newb =[b for b,bb in zip(bins,h) if subrange[0]<b<subrange[1]]
    newh =[bb for b,bb in zip(bins,h) if subrange[0]<b<subrange[1]]

    optimizedParams, pCov = opt.curve_fit(gaus, newb,newh, p0=[100, 1700, 100])
    plt.plot(newb,gaus(newb,*optimizedParams),'r--',label='fit', linewidth=2)
    # eventi sotto il picco
    C = optimizedParams[0]
    sigma = abs(optimizedParams[2])
    bw = bins[1]-bins[0]
    peak_events = C*sigma*np.sqrt(2*np.pi) / bw
    corrf = compute_corr_factor(nrun)
    rate = peak_events *  corrf / elapsed_time
    print('Run nr %d: events under the peak %.1f, correct. factor %.1f => rate %.f Hz' % (nrun, peak_events,corrf,rate))
    return rate

################
## MAIN PROGRAM
################

datadir  = './dati/'

# run list
runs = range(20272,20277)

rates=[]
for r in runs:
    rates.append(read_data_and_compute_rate(r))

Run nr 20272: events under the peak 17587.5, correct. factor 1.1 => rate 428 Hz
Run nr 20273: events under the peak 14271.8, correct. factor 1.0 => rate 141 Hz
Run nr 20274: events under the peak 11616.1, correct. factor 1.0 => rate 46 Hz
Run nr 20275: events under the peak 8611.0, correct. factor 1.0 => rate 16 Hz
Run nr 20276: events under the peak 5734.6, correct. factor 1.0 => rate 6 Hz

Usando i dati dei fit precedenti, fai il plot dei conteggi in funzione della distanza e fit con esponenziale

In [2]:
thick = [0,1,2,3,4] # cm

def expo(x, C, tau):
    return C*np.exp(-x*tau)


plt.plot(thick,rates,'bo', label='data')
oParams, pCov = opt.curve_fit(expo,thick,rates)
print(oParams)
plt.plot(np.linspace(0,4),expo(np.linspace(0,4),*oParams),'r--',label='fit', linewidth=2)
plt.legend()


## from mucal http://csrri.iit.edu/cgi-bin/mucal-form?name=Pb&ener=662000
expected=1.21584332

[428.09568958   1.11034418]
In [ ]: