Commit 92e29bf9 authored by Alfonso's avatar Alfonso
Browse files

dati mere gennaio

parent 28e986ff
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -102,4 +102,4 @@ def generate_fits(energy, theta_deg, phi_deg):
    except Exception as e:
        print("Errore durante la creazione del file .fits")

generate_fits(66, 0, 0)
 No newline at end of file
generate_fits(66, 40, 40)
 No newline at end of file
+45 −15
Original line number Diff line number Diff line
@@ -6,7 +6,7 @@ import astropy.io.fits as pyfits
from noisy import add_noise_X, add_noise_S


def write_fits_dep(energy, en_primary, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_id, x_pixel, y_pixel):  # ! è qui che decido quali informazioni sulla simulazione voglio nei file fits_dep
def write_fits_dep(event_id, energy, en_primary, theta, phi, type, en_dep_noisy, en_dep, name_id, pos_id, x_pixel, y_pixel):  # ! è qui che decido quali informazioni sulla simulazione voglio nei file fits_dep
    
    x1_pixel = (pos_id // 89) +1
    y1_pixel = (pos_id % 89) +1
@@ -19,13 +19,14 @@ def write_fits_dep(energy, en_primary, theta, phi, type, en_dep_noisy, en_dep, n
    
    # Extension
    tbhdu = pyfits.BinTableHDU.from_columns([
        #pyfits.Column(name='En_dep_Noise', format='1D', unit='keV', array = en_dep_noisy),
        pyfits.Column(name='EventID', format='1J', array = event_id),
        pyfits.Column(name='En_Primary', format='1D', unit='keV', array = en_primary),
        pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep),
        pyfits.Column(name='En_dep_Noise', format='1D', unit='keV', array = en_dep_noisy),
        #pyfits.Column(name=name_id, format='1I', array = pos_id),
        pyfits.Column(name='X_Y', format='2J', array = X_Y), 
        pyfits.Column(name='X_Y', format='2I', array = X_Y), 
        pyfits.Column(name='X_Pixel_cm', format='1D', unit='cm', array = x_pixel),
        pyfits.Column(name='Y_Pixel_cm', format='1D', unit='cm', array = y_pixel)
        #pyfits.Column(name='En_dep', format='1D', unit='keV', array = en_dep)
        ])
    
    tbhdu.header['EXTNAME'] = 'EVENTS'
@@ -66,35 +67,64 @@ def get_data(files):
    with pyfits.open(files) as hdul:
        data = hdul[1].data
        tot_events = len(data)
        
        pixel_id_col = data['Pixel_ID']
        pixel_id = pixel_id_col[pixel_id_col != -1000]
        
        scint_id_col = data['Scint_ID']
        scint_id = scint_id_col[scint_id_col != -1000]
        
        # mask_X = (data['Scint_ID'] == -1000)   # condizione              //TOdO questo è un altro modo di gestire i dati in modo piu compatto, non è concluso
        # eventi_X = data[mask_X]
        # mask_S = (data['Scint_ID'] != -1000)   # condizione
        # eventi_S = data[mask_S]
        # columns_for_X = ['EventID','En_Dep','Pixel_ID','X_Pixel','Y_Pixel','En_Primary']   #  qui scelgo le info
        # columns_for_S = ['EventID','En_Dep','Scint_ID','X_Pixel','Y_Pixel','En_Primary'] 
        # eventi_X_reduced = eventi_X[columns_for_X]
        # eventi_S_reduced = eventi_S[columns_for_S]
        
        event_id = data['EventID']
        X_event_id = event_id[scint_id_col == -1000] 
        S_event_id = event_id[scint_id_col != -1000]
        
        energy = data['En_Primary'][0]
        
        en_primary = data['En_Primary']
        X_en_primary = en_primary[scint_id_col == -1000] 
        S_en_primary = en_primary[scint_id_col != -1000]
        
        en_dep = data['En_dep']
        x_pixel = data['X_Pixel']
        y_pixel = data['Y_Pixel']
        X_en_dep = en_dep[scint_id_col == -1000]
        S_en_dep = en_dep[scint_id_col != -1000]
        eventi_X = (scint_id_col == -1000).sum()
        eventi_S = tot_events - eventi_X
        ratio_X = eventi_X/tot_events
        ratio_S = eventi_S/tot_events
        
        x_pixel = data['X_Pixel']
        X_x_pixel = x_pixel[scint_id_col == -1000]
        S_x_pixel = x_pixel[scint_id_col != -1000]
        
        y_pixel = data['Y_Pixel']
        X_y_pixel = y_pixel[scint_id_col == -1000]
        S_y_pixel = y_pixel[scint_id_col != -1000]
        
        num_eventi_X = (scint_id_col == -1000).sum()
        num_eventi_S = (scint_id_col != -1000).sum()
        ratio_X = num_eventi_X/tot_events
        ratio_S = num_eventi_S/tot_events
        if num_eventi_X + num_eventi_S != tot_events:
            raise ValueError("Errore: Eventi X + Eventi S DIVERSO da Eventi Totali!")

        X_en_dep_noisy = [add_noise_X(X_dep) for X_dep in X_en_dep] # list comprehension di Python, una struttura compatta che consente di applicare una funzione a ciascun elemento di un array o di una lista, creando direttamente un nuovo array con i risultati.
        S_en_dep_noisy = [add_noise_S(S_dep) for S_dep in S_en_dep]

    return energy, en_primary, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id, x_pixel, y_pixel
    return X_event_id, S_event_id, energy, X_en_primary, S_en_primary, ratio_X, ratio_S, X_x_pixel, S_x_pixel, X_y_pixel, S_y_pixel, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id

# scrive sia il file di testo che i file fits
def write_results(file_dat, results):

    with open(file_dat, "w") as f:
        for energy, en_primary, ratio_X, ratio_S, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id, x_pixel, y_pixel in results:
        for X_event_id, S_event_id, energy, X_en_primary, S_en_primary, ratio_X, ratio_S, X_x_pixel, S_x_pixel, X_y_pixel, S_y_pixel, theta, phi, X_en_dep_noisy, S_en_dep_noisy, X_en_dep, S_en_dep, pixel_id, scint_id in results:
            
            X_fits = write_fits_dep(energy, en_primary, theta, phi, 'X', X_en_dep_noisy, X_en_dep, 'Pixel_ID', pixel_id, x_pixel, y_pixel)
            S_fits = write_fits_dep(energy, en_primary, theta, phi, 'S', S_en_dep_noisy, S_en_dep, 'Scint_ID', scint_id, x_pixel, y_pixel)
            X_fits = write_fits_dep(X_event_id, energy, X_en_primary, theta, phi, 'X', X_en_dep_noisy, X_en_dep, 'Pixel_ID', pixel_id, X_x_pixel, X_y_pixel)
            S_fits = write_fits_dep(S_event_id, energy, S_en_primary, theta, phi, 'S', S_en_dep_noisy, S_en_dep, 'Scint_ID', scint_id, S_x_pixel, S_y_pixel)
            
            #f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\t{theta}\t{phi}\t{X_fits}\t{S_fits}\n")