Commit 66c74b6c authored by Alfonso's avatar Alfonso
Browse files

file .arf e plot

parent 2a10a1b5

python/aeFF.py

+49 −55
Original line number Diff line number Diff line 
# va bene fare la media per ogni combinazione theta phi di energia?

# dato che divido per gli eventi detected va bene cosi o devo dividere per il totale

# dei generati e moltiplicare per l'area della sorgente?

import os

import astropy.io.fits as fits

import pandas as pd

from multiprocessing import Pool

# come faccio a sapere se il mio codice è robusto



import numpy as np

from data import data_estraction

from plot import plot_area

from arf import arf_file



def main():

    # 1 Estract data from fits file

    file_dat = "areaEfficace.dat"

    data_estraction(file_dat)

    

    #2 Filter data

    data = np.loadtxt(file_dat)

    

    energy = data[:, 0]

    ratio_X = data[:, 1]

    ratio_S = data[:, 2]

    theta = data[:, 3]

    phi = data[:, 4]

    

    # Ottieni tutte le combinazioni uniche di theta e phi

    angle_combinations = np.unique(np.column_stack((theta, phi)), axis=0)

    

    # Inizializza due liste per tracciare i dati per conti X e conti S

    plot_data_X = []

    plot_data_S = []

    

    # Cicla attraverso tutte le combinazioni uniche di theta e phi

    for combination in angle_combinations:  # estrae dall'array due valori e gli da il nome

        theta_val = combination[0]

        phi_val = combination[1]



        # Maschera booleana per filtrare i dati per la combinazione corrente di theta e phi

        # confronta i dati del file con i dati dell'array creato

        mask = (theta == theta_val) & (phi == phi_val)

        energy_filtered = energy[mask]

        ratio_X_filtered = ratio_X[mask]

        ratio_S_filtered = ratio_S[mask]

        

        # Salva i dati per il plot

        plot_data_X.append((energy_filtered, ratio_X_filtered, theta_val, phi_val))

        plot_data_S.append((energy_filtered, ratio_S_filtered, theta_val, phi_val))

    

        # 3 Arf

        arf_file(energy_filtered, ratio_X_filtered, ratio_S_filtered, theta_val, phi_val)



    # 4 Plot

    plot_area(plot_data_X, 'X', file_name="XGIS_AeffX.png")

    plot_area(plot_data_S, 'S', file_name="XGIS_AeffS.png")



# Funzione che prende in input un singolo file FITS da processare.

def get_data(file_path):

    with fits.open(file_path) as hdul:

        data = hdul[1].data

        tot_events = len(data)

        scint_id_col = data['Scint_ID']

        energy = data['En_Primary'][0]

        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

        

        return energy, ratio_X, ratio_S



def write_raw_results(output, results):

    with open(output, "a") as f:

        for energy, ratio_X, ratio_S in results:

            f.write(f"{energy:.1f}\t{ratio_X:.3f}\t{ratio_S:.3f}\n")



# se vuoi i dati rough commenta questa funzione

#def write_results(output_file):

#    # Questa riga legge il file .dat e crea un DataFrame

#    data = pd.read_csv(output_file, delim_whitespace=True, header=None, names=['energy', 'ratio_X', 'ratio_S'])

#    # Calcola la media per ogni valore unico di energia

#    result = data.groupby('energy').mean().reset_index()

#    # Approssima i risultati alla terza cifra decimale

#    result = result.round(3)

#    # Scrive il risultato nel file di output

#    result.to_csv(output_file, sep='\t', index=False, header=False)



# Funzione che processa tutti i file .fits presenti nella folder

def process_folder_parallel(folder, output):

    # Crea una lista di file FITS trovati nella folder. Sorted per leggerli in ordine

    files = sorted([os.path.join(folder, file) for file in os.listdir(folder) if file.endswith(".fits")], key=os.path.getmtime)

    

    # Crea un pool di processi che permette di eseguire funzioni in parallelo.

    # Usa il metodo map() del pool per eseguire la funzione get_data() su ciascun file FITS in parallelo.

    # I risultati vengono raccolti in una lista chiamata results.

    with Pool() as pool:

        results = pool.map(get_data, files) # qui files è l'argomento che passi a get_data

    

    write_raw_results(output, results)

    

    #write_results(output_file)



if __name__ == "__main__":

    folder = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits"

    output = "areaEfficace.dat"

    main()



    

    # Inizializza il file di output

    open(output, "w").close()



    process_folder_parallel(folder, output)

    
    plot_area(output)
 
 No newline at end of file

python/arf.py

+53 −45
Original line number Diff line number Diff line
@@ -2,49 +2,57 @@ import astropy.io.fits as pyfits 
import os

import numpy as np



try:

def create_header(extension):

    extension.header['EXTNAME'] = 'SPECRESP'

    extension.header['TELESCOP'] = 'THESEUS'

    extension.header['INSTRUME'] = 'XGIS'

    extension.header['FILTER'] = 'NONE'

    extension.header['HDUCLASS'] = 'OGIP'

    extension.header['HDUCLAS1'] = 'RESPONSE'

    extension.header['HDUCLAS2'] = 'SPECRESP'

    extension.header['HDUVERS'] = '1.1.0'

    # obsolete (for old software)

    extension.header['ARFVERSN'] = '1992A'

    extension.header['HDUVERS1'] = '1.0.0'

    extension.header['HDUVERS2'] = '1.1.0'



def arf_file(energy_filt, ratio_X_filt, ratio_S_filt, theta, phi):

    try:

        # Creazione di "Null" primary array

        primary_hdu = pyfits.PrimaryHDU()



    # Extension

    data = np.loadtxt("areaEfficace.dat", comments='#')

    energy = data[:, 0]

    ratio_X = data[:, 1]

    ratio_S = data[:, 2]

    

    # Creazione di BinTableHDU

    bin_table_hdu = pyfits.BinTableHDU.from_columns([

        pyfits.Column(name='ENERG_LO', format='1E', unit='keV', array = energy[:-1]), # [inizio:fine:passo]

        pyfits.Column(name='ENERG_HI', format='1E', unit='keV', array = energy[1:]),  # 0 primo, 1 secondo, -1 ultimo escluso

        pyfits.Column(name='SPECRESP', format='1E', unit='cm**2', array = ratio_X),

        # Creazione di BinTableHDU per EVENTI X

        bin_tableX_hdu = pyfits.BinTableHDU.from_columns([

            pyfits.Column(name='ENERG_LO', format='1E', unit='keV', array=energy_filt[:-1]),

            pyfits.Column(name='ENERG_HI', format='1E', unit='keV', array=energy_filt[1:]),

            pyfits.Column(name='SPECRESP', format='1E', unit='cm**2', array=ratio_X_filt),

        ])

        create_header(bin_tableX_hdu)  # Aggiunta header comune per gli eventi X



    # Modifica dell'header. Mandatory keywords

    bin_table_hdu.header['EXTNAME'] = 'SPECRESP'

    bin_table_hdu.header['TELESCOP'] = 'THESEUS'

    bin_table_hdu.header['INSTRUME'] = 'XGIS'

    bin_table_hdu.header['FILTER'] = 'NONE'

    bin_table_hdu.header['HDUCLASS'] = 'OGIP'

    bin_table_hdu.header['HDUCLAS1'] = 'RESPONSE'

    bin_table_hdu.header['HDUCLAS2'] = 'SPECRESP'

    bin_table_hdu.header['HDUVERS'] = '1.1.0'

    # obsolete (for old software)

    bin_table_hdu.header['ARFVERSN'] = '1992A'

    bin_table_hdu.header['HDUVERS1'] = '1.0.0'

    bin_table_hdu.header['HDUVERS2'] = '1.1.0'  

    # optional

    # bin_table_hdu.header['PHAFILE'] = 

        # Creazione di BinTableHDU per EVENTI S

        bin_tableS_hdu = pyfits.BinTableHDU.from_columns([

            pyfits.Column(name='ENERG_LO', format='1E', unit='keV', array=energy_filt[:-1]),

            pyfits.Column(name='ENERG_HI', format='1E', unit='keV', array=energy_filt[1:]),

            pyfits.Column(name='SPECRESP', format='1E', unit='cm**2', array=ratio_S_filt),

        ])

        create_header(bin_tableS_hdu)  # Aggiunta header comune per gli eventi S



        # Creazione di HDUList e scrittura del file FITS

    hdul = pyfits.HDUList([primary_hdu , bin_table_hdu])

        hdulX = pyfits.HDUList([primary_hdu, bin_tableX_hdu])

        hdulS = pyfits.HDUList([primary_hdu, bin_tableS_hdu])



        output_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/arf"

        

    output_dir = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python"

    #output_file = os.path.join(output_dir, input_file.split(".")[0] + ".arf")

    output_file = os.path.join(output_dir, "areaEfficace.arf")

        if not os.path.exists(output_dir):

            os.makedirs(output_dir)

            

    hdul.writeto(output_file, overwrite=True)

        output_X = os.path.join(output_dir, f"XaeFF_{theta}_{phi}.arf")

        output_S = os.path.join(output_dir, f"SaeFF_{theta}_{phi}.arf")



    print("File .arf creato correttamente")

        hdulX.writeto(output_X, overwrite=True)

        hdulS.writeto(output_S, overwrite=True)



        print("File .arf creati correttamente.")

    except Exception as e:

    print("Errore durante la creazione del file .arf")
 
 No newline at end of file
 
        print(f"Errore durante la creazione del file .arf: {e}")

python/plot.py

+15 −24
Original line number Diff line number Diff line
@@ -2,40 +2,31 @@ import numpy as np 
import matplotlib.pyplot as plt

import os



def plot_area(output):

    # Legge il file ignorando le righe che iniziano con #

    data = np.loadtxt(output, comments='#')

def plot_area(plot_data, type, file_name):



    # Separa le due colonne

    energy = data[:, 0]  # Prima colonna: Energy [keV]

    ratio_X = data[:, 1]

    ratio_S = data[:, 2]

    for data in plot_data:

        energy, ratio, theta_val, phi_val = data

        label = f'Theta={theta_val}, Phi={phi_val}'

        

    # Crea il grafico

    plt.figure(figsize=(8, 6))

    plt.plot(energy, ratio_X, marker='o', linestyle='-', color='b', label='sdd')

    plt.plot(energy, ratio_S, marker='x', linestyle='--', color='r', label='scint')

        plt.plot(energy, ratio, marker='o', linestyle='-', label=label)



    # Aggiungi etichette e titolo

    plt.xlabel('Energy [keV]')

    plt.ylabel('Aeff')

    plt.title('XGIS')

    plt.ylabel(f'Aeff_{type}')

    plt.title(f'XGIS - Eventi_{type}')

    plt.grid(True)



    # Mostra la legenda

    plt.legend()

    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')



    # Salva il grafico

    # Specifica il percorso completo dove salvare il grafico

    directory = '/home/alfonso/Scrivania/THESEUS/xgis_m7-main/python/figure'  

    if not os.path.exists(directory):

        os.makedirs(directory)



    file_path = os.path.join(directory, output.split(".")[0] + ".png")  #.pdf



    # Salva il grafico

    plt.savefig(file_path)



    # Mostra il grafico

    plt.show()

    file_path = os.path.join(directory, file_name)

    plt.savefig(file_path, bbox_inches='tight')

    

    # Pulire la figura per i plot successivi

    plt.clf() 

    #plt.show()
 
 No newline at end of file