From 66c74b6c4c99591cde42b30605a7ac1857c501f3 Mon Sep 17 00:00:00 2001
From: Alfonso <alo.99@hotmail.it>
Date: Fri, 4 Oct 2024 11:15:54 +0200
Subject: [PATCH] file .arf e plot

---
 python/aeFF.py | 104 +++++++++++++++++++++++--------------------------
 python/arf.py  |  98 +++++++++++++++++++++++++---------------------
 python/plot.py |  39 +++++++------------
 3 files changed, 117 insertions(+), 124 deletions(-)

diff --git a/python/aeFF.py b/python/aeFF.py
index ec8d17f..38d0a97 100644
--- a/python/aeFF.py
+++ b/python/aeFF.py
@@ -1,64 +1,58 @@
-# 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
-from plot import plot_area
-
-# 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")
+# come faccio a sapere se il mio codice è robusto
 
-# 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)
+import numpy as np
+from data import data_estraction
+from plot import plot_area
+from arf import arf_file
 
-# 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)
+def main():
+    # 1 Estract data from fits file
+    file_dat = "areaEfficace.dat"
+    data_estraction(file_dat)
+    
+    #2 Filter data
+    data = np.loadtxt(file_dat)
     
-    # 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
+    energy = data[:, 0]
+    ratio_X = data[:, 1]
+    ratio_S = data[:, 2]
+    theta = data[:, 3]
+    phi = data[:, 4]
     
-    write_raw_results(output, results)
+    # Ottieni tutte le combinazioni uniche di theta e phi
+    angle_combinations = np.unique(np.column_stack((theta, phi)), axis=0)
     
-    #write_results(output_file)
+    # 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]
 
-if __name__ == "__main__":
-    folder = "/home/alfonso/Scrivania/THESEUS/xgis_m7-main/fits"
-    output = "areaEfficace.dat"
+        # 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))
     
-    # Inizializza il file di output
-    open(output, "w").close()
+        # 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")
+
+
+if __name__ == "__main__":
+    main()
+
     
-    process_folder_parallel(folder, output)
+
     
-    plot_area(output)
\ No newline at end of file
diff --git a/python/arf.py b/python/arf.py
index ec6749b..69abc57 100644
--- a/python/arf.py
+++ b/python/arf.py
@@ -2,49 +2,57 @@ import astropy.io.fits as pyfits
 import os
 import numpy as np
 
-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),
-    ])
-    
-    # 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'
+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)
-    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 HDUList e scrittura del file FITS
-    hdul = pyfits.HDUList([primary_hdu , bin_table_hdu])
-    
-    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")
-    
-    hdul.writeto(output_file, overwrite=True)
-
-    print("File .arf creato correttamente")
-except Exception as e:
-    print("Errore durante la creazione del file .arf")
\ No newline at end of file
+    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()
+
+        # 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
+
+        # 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
+        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"
+        
+        if not os.path.exists(output_dir):
+            os.makedirs(output_dir)
+            
+        output_X = os.path.join(output_dir, f"XaeFF_{theta}_{phi}.arf")
+        output_S = os.path.join(output_dir, f"SaeFF_{theta}_{phi}.arf")
+
+        hdulX.writeto(output_X, overwrite=True)
+        hdulS.writeto(output_S, overwrite=True)
+
+        print("File .arf creati correttamente.")
+    except Exception as e:
+        print(f"Errore durante la creazione del file .arf: {e}")
+
diff --git a/python/plot.py b/python/plot.py
index 0030fa3..e3e9b16 100644
--- a/python/plot.py
+++ b/python/plot.py
@@ -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]
-
-    # 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')
+    for data in plot_data:
+        energy, ratio, theta_val, phi_val = data
+        label = f'Theta={theta_val}, Phi={phi_val}'
+        
+        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
-- 
GitLab