Commit 75356ee9 authored by Patrizia Barria's avatar Patrizia Barria
Browse files

First upload test

parent 8cb643c9
Loading
Loading
Loading
Loading
+473 −0
Original line number Original line Diff line number Diff line
import seaborn as sns
import os
import pandas as pd
import warnings
import ast
import pickle
#---- Paralelization ----
# import os
# os.environ["MODIN_ENGINE"] = "ray"
# import ray
# ray.init()
# import modin.pandas as pd
#++++++++++++++++++++++++++
# import os
# os.environ["RAY_ENABLE_MAC_LARGE_OBJECT_STORE"] = "1"
# import ray
# # print(ray.__version__)
# ray.init(object_store_memory=2 * 1024 * 1024 * 1024)  # 2GB
# import modin.pandas as pd
# import modin
# # print(modin.__version__)
#++++++++++++++++++++++++
from array import array
import sys, os, argparse
import math as mt
import numpy as np
from numpy import linspace
import re
from collections import defaultdict
import csv
import json
import itertools
from prettytable import PrettyTable
import argparse
import glob
from adjustText import adjust_text
from IPython.display import display, HTML
import pywt
import math
import peakutils
import itertools
import statsmodels.api as sm
from IPython.display import display, clear_output
from iminuit import Minuit
from iminuit.cost import LeastSquares
from skimage.filters import threshold_otsu
from sklearn.mixture import GaussianMixture
from sklearn.linear_model import LinearRegression, RANSACRegressor
#++++ ASTROPY +++++++++++
from astropy.io import fits
from astropy.table import Table
from astropy.visualization import hist
from astropy import stats
from astropy.stats import scott_bin_width
#++++ MATPLOTLIB +++++++++++
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as pp
import matplotlib.patheffects as path_effects
import matplotlib.cm as cm
from matplotlib import animation, rc
from pathlib import Path
from matplotlib import rc, pyplot as plt
from matplotlib.ticker import PercentFormatter
from matplotlib.patches import Rectangle
from matplotlib import rc
from matplotlib.gridspec import GridSpec
from matplotlib.pylab import *
from matplotlib.animation import ArtistAnimation
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
import matplotlib.animation as animation
from matplotlib.ticker import ScalarFormatter
from matplotlib.animation import PillowWriter
from matplotlib import cm  # Import the cm module for colormaps
from matplotlib.animation import FuncAnimation, FFMpegWriter
from matplotlib.colors import Normalize
from matplotlib.colors import LogNorm
from matplotlib.cm import get_cmap
import matplotlib.ticker as ticker
from matplotlib.ticker import NullFormatter, MaxNLocator
from matplotlib.colors import ListedColormap
#%matplotlib inline
#%time
#++++ SCIPY ++++++++++++++++
import scipy.optimize as opt
import scipy.stats as stat
from scipy.stats import iqr
from scipy.interpolate import interp1d
from scipy.stats import norm
from scipy.stats import poisson
from scipy.stats import expon
from scipy.signal import savgol_filter
from scipy.stats import *
from scipy import interpolate
from scipy.interpolate import interp2d
import statsmodels.api as sm
import scipy.signal as ss
from scipy.optimize import minimize_scalar
from scipy.ndimage.filters import maximum_filter
from scipy.ndimage.morphology import generate_binary_structure, binary_erosion
from scipy.signal import find_peaks, peak_widths
from scipy.ndimage.interpolation import shift
from scipy.fft import fft
from scipy.optimize import curve_fit, minimize
from scipy.signal import spectrogram
from scipy.signal import deconvolve
from scipy.signal import wiener
from scipy.signal import *
from scipy import signal
from scipy.special import voigt_profile
from scipy.stats import zscore
from scipy.stats import gaussian_kde
#++++ DATETIME +++++++++++++
import datetime
from datetime import datetime, timedelta
import time
#++++ MPL TOOLKITS +++++++++
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
from mpl_toolkits.axes_grid1 import make_axes_locatable



#+++++++++++++++++++++++++++++++++++++++++++++++++
#            OUTPUT DIRECTORY
#+++++++++++++++++++++++++++++++++++++++++++++++++

# output_directory3 = '/Users/merlino/Desktop/INAF-IAPS/20231025_camera_climatica_cs137/TEST/Plots/' #Cs137
# output_directory4 = '/Users/merlino/Desktop/INAF-IAPS/20231026_camera_climatica_ba133/TEST/Plots/' #Ba133
# output_directory5 = '/Users/merlino/Desktop/INAF-IAPS/20240227_camera_climatica_lyso/TEST/Plots/' #LYSO
# output_directory6 = '/Users/merlino/Desktop/INAF-IAPS/20240228_camera_climatica_fondo/TEST/Plots/' #BKG
output_directory3 = '/Users/merlino/Desktop/INAF-IAPS/20240618_camera_climatica/co57/TEST/Plots/' #Am241
# output_directory3 = '/Users/merlino/Desktop/INAF-IAPS/20240618_camera_climatica/fondo/TEST/Plots/' #BKG2
#output_directory3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/volo/TEST/Plots/'  # Timmins Flight
#output_directory3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/post_volo/TEST/Plots/'  # Timmins Post Flight


#++++++Commom Output Directory++++++++++++++++++
# output_directory = '/Users/merlino/Desktop/INAF-IAPS/GRASSII_2023-2024-Thermal_Cycles/Ba133_Cs137_LYSO/Event_Matching/Cs137/'
#output_directory = '/Users/merlino/Desktop/INAF-IAPS/GRASSII_2023-2024-Thermal_Cycles/Am241-Ba133_Cs137_Co57_Na22/Co57/Event_Matching/Anger_Data/'
output_directory = '/Users/merlino/Desktop/INAF-IAPS/GRASSII_2023-2024-Thermal_Cycles/Co57/Event_Matching/Anger_Data/'
#output_directory = '/Users/merlino/Desktop/INAF-IAPS/GRASSII_2023-2024-Thermal_Cycles/Timmins_Fight/Event_Matching/Anger_Data/'
#output_directory = '/Users/merlino/Desktop/INAF-IAPS/GRASSII_2023-2024-Thermal_Cycles/Timmins_Post-Fight/Event_Matching/Anger_Data/'

#+++++ Create Output Directory if doesn't exist
os.makedirs(output_directory, exist_ok=True)
os.makedirs(output_directory3, exist_ok=True)
# os.makedirs(output_directory4, exist_ok=True)
# os.makedirs(output_directory5, exist_ok=True)
# os.makedirs(output_directory6, exist_ok=True)



#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#             FILES TO PROCESS
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++

#---> Directories where CSV files are located
# columns_to_exclude = ['ec'] # 'ec' has to be excluded for Tammins files
columns_to_exclude = [''] # 'ec' has to be excluded for Tammins files

#------------------> SIGNALS

#---> Cs137
# directory_path3 = '/Users/merlino/Desktop/INAF-IAPS/20231025_camera_climatica_cs137/adc/'
#---> Ba133
# directory_path4 = '/Users/merlino/Desktop/INAF-IAPS/20231026_camera_climatica_ba133/adc/'
#---> LYSO
# directory_path5 = '/Users/merlino/Desktop/INAF-IAPS/20240227_camera_climatica_lyso/adc/'
#---> Am241
directory_path3 = '/Users/merlino/Desktop/INAF-IAPS/20240618_camera_climatica/co57/adc/'
#---> Timmins Flight
#directory_path3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/volo/adc/'
#directory_path3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/post_volo/adc/'

#------------------> BACKGROUNDS
# directory_path6 = '/Users/merlino/Desktop/INAF-IAPS/20240228_camera_climatica_fondo/adc/'


#------------------> REGISTERS

#---> Cs137
# temperature_data_file3 = '/Users/merlino/Desktop/INAF-IAPS/20231025_camera_climatica_cs137/TEST/combined_registers.csv'
#---> Ba133
# temperature_data_file4 = '/Users/merlino/Desktop/INAF-IAPS/20231026_camera_climatica_ba133/TEST/combined_registers.csv'
#---> LYSO
# temperature_data_file5 = '/Users/merlino/Desktop/INAF-IAPS/20240227_camera_climatica_lyso/TEST/combined_registers.csv'
#---> BACKGROUNDS
# temperature_data_file6 = '/Users/merlino/Desktop/INAF-IAPS/20240228_camera_climatica_fondo/TEST/combined_registers.csv'
#---> Am241
temperature_data_file3 = '/Users/merlino/Desktop/INAF-IAPS/20240618_camera_climatica/co57/TEST/combined_registers.csv'
#---> Timmins Flight
#temperature_data_file3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/volo/TEST/combined_registers.csv'
#temperature_data_file3 = '/Users/merlino/Desktop/INAF-IAPS/Flight_campaign_TIMMINS/post_volo/TEST/combined_registers.csv'


#++++++++++++++++++++++++++++++
#    ARRAY DIMENSIONS
#++++++++++++++++++++++++++++++
# height = 50.44 mm
# height = (6.13*8+(0.2*7) mm
# active area = 6.07 mm
# pitch = 6.13 mm +/- 0.05
# gap = 0.2 mm
#++++++++++++++++++++++++++++++


# Define constants
CRYSTAL_WIDTH_MM = 49.04  # Crystal dimensions should be 50.44 
CRYSTAL_HEIGHT_MM = 49.04
#CRYSTAL_WIDTH_MM = 50.44  # Crystal dimensions should be 50.44 
#CRYSTAL_HEIGHT_MM = 50.44

SIPM_ARRAY_WIDTH_MM = 6.07 * 8  # SiPM array dimensions
SIPM_ARRAY_HEIGHT_MM = 6.07 * 8

# Conversion factors
X_CONVERSION_FACTOR = CRYSTAL_WIDTH_MM / 8
Y_CONVERSION_FACTOR = CRYSTAL_HEIGHT_MM / 8

# Fractions for X and Y
X_MINUS_FRACTIONS = [0.1250, 0.2483, 0.3731, 0.5000, 0.6250, 0.7500, 0.8876, 1.0000]
X_PLUS_FRACTIONS = [1.1250 - fraction for fraction in X_MINUS_FRACTIONS]
Y_MINUS_FRACTIONS = [1.0000, 0.8876, 0.7500, 0.6250, 0.5000, 0.3731, 0.2483, 0.1250]
Y_PLUS_FRACTIONS = [1.1250 - fraction for fraction in Y_MINUS_FRACTIONS]

def process_file(file_path, output_directory3):


#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#             ADDING EXECUTION TIME 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    start = time.time()
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++

    # Read DataFrame
    df = pd.read_csv(file_path)

        
    filename = os.path.basename(file_path)

    # Extract temperature value
    temperature = df['Temperature'].iloc[0]
    print(f"Processing file: {filename} @T={temperature}°C: {len(df)} entries")

    #------- Here we save only Anger coordinates ---------------
#     # Initialize an empty dictionary to store channel pairs
#     channel_pairs = define_channel_pairs()

#     # Compute the sum of signals for each pixel
#     pixel_signal_df = compute_pixel_signals(df, channel_pairs)

#     # Process pixel signals to compute various metrics
#     df_single = process_pixel_signals(df, pixel_signal_df)

#     # Save modified DataFrame to CSV
#     save_to_csv(df_single, file_path, output_directory)

#     # Print processing summary
#     print_processing_summary(df, df_single)
    
    #------- Here we append Anger coordinates to original dataframe ---------------
    
    # Initialize an empty dictionary to store channel pairs
    channel_pairs = define_channel_pairs()

    # Compute the sum of signals for each pixel
    pixel_signal_df = compute_pixel_signals(df, channel_pairs)

    # Process pixel signals to compute various metrics
#     df_single = process_pixel_signals(df, pixel_signal_df)
    df_single = process_pixel_signals(df, pixel_signal_df, channel_pairs)


    # Append computed metrics to the original DataFrame
    df = pd.concat([df, df_single], axis=1)

    # Save modified DataFrame to CSV
    save_to_csv(df, file_path, output_directory)

    # Print processing summary
    print_processing_summary(df, df_single)
    

#++++++++++++++++++++++++++++++++++++++++++
#      PRINT OUT TIME TO PROCESS SCRIPT
#++++++++++++++++++++++++++++++++++++++++++

    # Capture end time and calculate the time taken
    end = time.time()
    time_taken = end - start
    print(f"Time taken to process {filename}: {time_taken}\n")

    # Convert the time taken into hours, minutes, and seconds
    hours, minutes, seconds = convert_seconds_to_hms(time_taken)

    # Print the result
    print(f"Time taken to process {filename}: {hours} hours, {minutes} minutes, and {seconds:.2f} seconds")

#++++++++++++++++++++++++++++++++++++++++++

def convert_seconds_to_hms(seconds):
    # Convert total seconds to hours, minutes, and remaining seconds
    hours = seconds // 3600
    minutes = (seconds % 3600) // 60
    remaining_seconds = seconds % 60
    return int(hours), int(minutes), remaining_seconds    


def define_channel_pairs():
    channel_pairs = {}
    for pixel_num in range(1, 65):
        pixel_name = f'pixel{pixel_num}'
        channel1 = f'ch{((pixel_num - 1) % 8) + 1}'
        channel2 = f'ch{8 + ((pixel_num - 1) // 8) + 1}'
        channel_pairs[pixel_name] = [channel1, channel2]
    return channel_pairs

def compute_pixel_signals(df, channel_pairs):
    pixel_signal_df = pd.DataFrame()
    for pixel, channels in channel_pairs.items():
        pixel_signal_df[pixel] = df[channels[0]] + df[channels[1]]
    return pixel_signal_df

# def process_pixel_signals(df, pixel_signal_df):
def process_pixel_signals(df, pixel_signal_df, channel_pairs):    
    # Initialize lists for storing data
    single_event = []
    sig_cross_event = []
    x_physical_list = []
    y_physical_list = []
    x_column_list = []
    y_row_list = []
    row_channel_number_list = []
    col_channel_number_list = []
    row_channel_value_list = []
    col_channel_value_list = []

    # Iterate over rows and compute various metrics
    for idx, row in df.iterrows():
        matrix_df = pd.DataFrame(index=[f'ch{i + 1}' for i in range(8)], columns=[f'ch{i + 9}' for i in range(8)])

        # Calculate matrix values using pixel_signal_df
        for pixel, channels in channel_pairs.items():
            matrix_df.loc[channels[0], channels[1]] = pixel_signal_df.loc[idx, pixel]

        # Calculate maximum value and its coordinates
        max_value = matrix_df.values.max()
        max_coords = np.unravel_index(matrix_df.values.argmax(), matrix_df.shape)

        # Calculate signal cross

        sig_cross = matrix_df.iloc[max_coords[0], :].sum() + matrix_df.iloc[:, max_coords[1]].sum() - max_value
        sig_cross_event.append(sig_cross)

        # Calculate physical coordinates
        x_physical = max_coords[1] * X_CONVERSION_FACTOR
        y_physical = max_coords[0] * Y_CONVERSION_FACTOR

        # Calculate row and column channel numbers and values
        row_channel_number = max_coords[0] + 1
        col_channel_number = max_coords[1] + 9
        row_channel_value = row[f'ch{row_channel_number}']
        col_channel_value = row[f'ch{col_channel_number}']
        
        # Convert row_channel_value and col_channel_value to integers
        row_channel_value = int(row_channel_value)
        col_channel_value = int(col_channel_value)

        
#         print(row_channel_number)
#         print(col_channel_number)
#         print(row_channel_value)
#         print(col_channel_value)
#         print(max_value)
#         print(sig_cross)

        # Calculate fractions for X and Y
        x_minus_fraction = X_MINUS_FRACTIONS[max_coords[0]]
        x_plus_fraction = X_PLUS_FRACTIONS[max_coords[0]]
        y_minus_fraction = Y_MINUS_FRACTIONS[max_coords[1]]
        y_plus_fraction = Y_PLUS_FRACTIONS[max_coords[1]]

        # Calculate X column and Y row
        x_column = round((x_plus_fraction - x_minus_fraction) / (x_plus_fraction + x_minus_fraction), 3)
        y_row = round((y_plus_fraction - y_minus_fraction) / (y_plus_fraction + y_minus_fraction), 3)

        # Append values to lists
        single_event.append(max_value)
        x_physical_list.append(x_physical)
        y_physical_list.append(y_physical)
        x_column_list.append(x_column)
        y_row_list.append(y_row)
        row_channel_number_list.append(row_channel_number)
        col_channel_number_list.append(col_channel_number)
        row_channel_value_list.append(row_channel_value)
        col_channel_value_list.append(col_channel_value)

    # Create DataFrame for processed metrics
    df_single = pd.DataFrame({
        'single_event': single_event,
        'sig_cross_event': sig_cross_event,
        'x_physical_list': x_physical_list,
        'y_physical_list': y_physical_list,
        'x_column': x_column_list,
        'y_row': y_row_list,
        'row_channel_number': row_channel_number_list,
        'col_channel_number': col_channel_number_list,
        'row_channel_value': row_channel_value_list,
        'col_channel_value': col_channel_value_list
    })

    return df_single

def save_to_csv(df_single, file_path, output_directory):
    # Define output file path
    output_file = os.path.join(output_directory, os.path.basename(file_path))

    # Save modified DataFrame to CSV
    try:
        df_single.to_csv(output_file, index=False)
        print(f"#---> Saved {output_file} successfully")
    except Exception as e:
        print(f"#---> Error saving {output_file}: {e}")

def print_processing_summary(df, df_single):
    # Print processing summary
    print(f"Finished processing DataFrame: {len(df)} rows")
    print(100 * '-')
    print('++++++ SIG w/signal and w/single event ++++++')
    print(100 * '-')
    print(df)
#     print(df_single)
    print(100 * '-')



def main():

#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#             ADDING EXECUTION TIME 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#    startTime = datetime.now()
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++


#     output_directory = "output_directory"
#     output_directory3 = "output_directory3"

    for filename in os.listdir(output_directory3):
        if filename.endswith('_data.csv'):
            file_path = os.path.join(output_directory3, filename)
#            process_file(file_path, output_directory)

       	    # Check if the processed file already exists in the output directory
            processed_file_path = os.path.join(output_directory, filename)
            if os.path.exists(processed_file_path):
               print(f"Skipping file {filename}: already processed")
               continue
        
       	    process_file(file_path, output_directory)

#++++++++++++++++++++++++++++++++++++++++++
#      PRINT OUT TIME TO PROCESS SCRIPT
#++++++++++++++++++++++++++++++++++++++++++
  # Print the time taken to process all files
#    print(f"Total processing time: {datetime.now() - startTime}")
#++++++++++++++++++++++++++++++++++++++++++

if __name__ == "__main__":
    main()