diff --git a/Electron_detector_eff_analyzer.py b/Electron_detector_eff_analyzer.py new file mode 100644 index 0000000..d8f112b --- /dev/null +++ b/Electron_detector_eff_analyzer.py @@ -0,0 +1,110 @@ +import numpy as np +import pandas as pd +import matplotlib.pylab as plt +from matplotlib.colors import ListedColormap +import seaborn as sns +import os +plt.style.use("IceCube") +import matplotlib as mpl +import numpy as np + +class DECOLeptonAnalyzer(): + r'''This class is for making plots like the ones in the + notebook that read in a bunch of simulated files + and make analysis level plots''' + + def __init__(self, pid, phi, thickness): + self.pid = pid + + self.thichness = thickness + + self.phi = phi + + + def read_hit_file(self, filename): + f = open(filename, 'r') + + xhits, yhits, charge = [], [], [] + + # skip first 2 lines + f.readline() + f.readline() + + x, y, c = [], [], [] + + flag = 1 + while 1: + + temp = f.readline().split() + + if len(temp) < 1 or temp[0] == '#': + break + + if temp[0] == '===': + continue + if temp[0] == '---': + if flag == 1: + flag = 0 + else: + xhits.append(x) + yhits.append(y) + charge.append(c) + x, y, c = [], [], [] + + else: + x.append(float(temp[1][:-1])) + y.append(float(temp[2][:-1])) + c.append(float(temp[3][:-1])) + + if len(x) > 0: + xhits.append(x) + yhits.append(y) + charge.append(c) + + return xhits, yhits, charge + + + + + def get_E_deposited(self, en, ang, posz): + x, y, c = self.read_hit_file( + "./output/{}/{}_posz_{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, round(posz, 3), en, float(ang), + float(self.phi), self.thichness)) + + dE_list = [] + + for i in range(len(c)): + dE_list.append(np.array(c[i]).sum() * 3.62) + + dE_list = np.array(dE_list) + + return dE_list.mean(), dE_list.std() + + + +energy = "10keV" + +theta = '90' +phi = '0' +particle_type = 'e-' + +dep_thickness = 26.3 + +z_pos = np.linspace(-1 * dep_thickness/2, dep_thickness/2, 100) + +a = DECOLeptonAnalyzer(particle_type, phi, dep_thickness) + +E_list = [] +err_list = [] + +for z in z_pos: + mean_val, std_val = a.get_E_deposited(energy, theta, z) + + E_list.append(mean_val) + err_list.append(std_val) + +f = plt.figure() +plt.errorbar(z_pos, E_list, yerr=err_list) +plt.xlabel("initial z position in um") +plt.ylabel("calibrated E in eV") +plt.show() \ No newline at end of file diff --git a/Electron_detector_eff_test.py b/Electron_detector_eff_test.py new file mode 100644 index 0000000..10acf03 --- /dev/null +++ b/Electron_detector_eff_test.py @@ -0,0 +1,164 @@ +r'''Class that handles various inputs and systematic parameters + and simulates muons''' +import os, sys +import subprocess +from pipes import quote +from math import * +import numpy as np + +class DECOMuonSimulator(): + r'''Simulator class for muons that uses allpix + and GEANT''' + + def __init__(self, pid, energy, pos, theta, **kwargs): + self.pid = pid #Particle id, ie mu+, e- + self.pos = pos + self.energy = energy + self.theta = float(theta) + self.phi = float(kwargs.pop('phi', 0.)) + self.depletion_thickness = float(kwargs.pop('depletion_thickness', 26.3)) + #todo: NOTE: in unit of um !!!!!!! + # Set other possible systematics with kwargs, including + # electric fields, pixel size, etc. + + self.path_to_root = os.getenv('DECO_ROOT_PATH') + self.path_to_geant = os.getenv('DECO_GEANT_PATH') + self.path_to_allpix = os.getenv('DECO_ALLPIX_PATH') + self.base_command = self.path_to_allpix + ' -c ./htc_wildfire/source_measurement.conf -o' + + def write_source_file(self, n_events, want_plot): + + with open('./htc_wildfire/source_measurement_replace.conf', 'r') as f: + data = f.readlines() + data[3] = data[3].format(n_events) + + data[18] = data[18].format(str(self.pos[0]) + " " + str(self.pos[1]) + " " + str(self.pos[2]) + "um") + + theta = pi - radians(self.theta) + phi = radians(self.phi) + dirx = 1*sin(theta)*cos(phi) + diry = 1*sin(theta)*sin(phi) + dirz = 1*cos(theta) + + data[20] = data[20].format(str(dirx) + " " + str(diry) + " " + str(dirz)) + + if want_plot == 'true': + data[38] = data[38].format('true' + "\n" + "output_linegraphs = true \n output_plots_step = 100ps \n output_plots_align_pixels = true \n output_plots_use_pixel_units = true") + else: + data[38] = data[38].format("false") + + + with open('./htc_wildfire/source_measurement.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + """ + def write_detector_file(self): + # USE THE PARAMETERS TO REWRITE DETECTOR FILE + with open('./htc_wildfire/detector_replace.conf', 'r') as f: + data = f.readlines() + data[-2] = data[-2].format(self.phi, self.theta, 0) + with open('./htc_wildfire/detector.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + with open('./htc_wildfire/htc_wildfire_shielded_replace.conf', 'r') as f: + data = f.readlines() + data[4] = data[4].format(self.depletion_thickness) + with open('./htc_wildfire/htc_wildfire_shielded.conf', 'w') as wf: + wf.writelines(data) + wf.close() + """ + + def set_output_file_name(self): + # write unique file name depending on parameters + if not os.path.exists("./output"): + os.system("mkdir ./output") + if not os.path.exists("./output/" + str(self.pid)): + os.system("mkdir ./output/" + str(self.pid)) + + self.outfile = str(self.pid) + "/" + "{}_posz_{}_theta_{:.1f}_phi_{:.1f}_thickiness_{:.1f}_highstats.txt".format(round(self.pos[2], 3), self.energy, self.theta, self.phi, self.depletion_thickness) + pass + + + def get_output_file_name(self): + try: + return self.outfile + except: + self.set_output_file_name() + return self.outfile + + + def run_simulation(self, n_events, want_charge_plot='false'): + + self.source_local_env() + # Run the allpix simulation + + output_file = self.get_output_file_name() + self.write_source_file(n_events, want_charge_plot) + + my_command = self.base_command[:] + + my_command += ' DepositionGeant4.particle_type="{}"'.format(self.pid) + my_command += ' -o DepositionGeant4.source_energy="{}"'.format(self.energy) + my_command += ' -o TextWriter.file_name="' + output_file + '"' + + """check if simulated using file name""" + if self.check_if_simulated(output_file) is True: + print("has been simulated") + return + + subprocess.call(my_command, shell=True) + + return + + + def check_if_simulated(self, filename): + # Check to see if simulation has already been run for + # this set of parameters + + if not os.path.exists("./output/" + filename): + return False + + return True + + + + def source_local_env(self): + geant = self.path_to_geant + "/bin/geant4.sh" + root = self.path_to_root + "/bin/thisroot.sh" + + pass + # Run source scripts for ROOT and GEANT if the + # simulation doesn't work + + +dep_thickness = 26.3 #um +pos_list = [] + +z_pos = np.linspace(-1 * dep_thickness/2, dep_thickness/2, 100) + +for i in range(len(z_pos)): + pos_list.append([0, 0, z_pos[i]]) + +#angles = ['0', '15', '30', '45', '60', '75'] + +#phis = ['0', '15', '30', '45', '60', '75', '90'] + +energy = ['10keV'] +angles = ['90'] +phis = ['0'] +particle_type = 'e-' + + +want_charge_plot = "false" + +for ene in energy: + for ang in angles: + for azi in phis: + for pos in pos_list: + a = DECOMuonSimulator(particle_type, ene, pos, ang, phi=azi, depletion_thickness=str(dep_thickness)) + + a.run_simulation(100, want_charge_plot) + + diff --git a/Electron_energy_caliberate.py b/Electron_energy_caliberate.py new file mode 100644 index 0000000..3a33510 --- /dev/null +++ b/Electron_energy_caliberate.py @@ -0,0 +1,166 @@ +r'''Class that handles various inputs and systematic parameters + and simulates muons''' +import os, sys +import subprocess +from pipes import quote +from math import * +import numpy as np + +class DECOMuonSimulator(): + r'''Simulator class for muons that uses allpix + and GEANT''' + + def __init__(self, pid, energy, pos, theta, **kwargs): + self.pid = pid #Particle id, ie mu+, e- + self.pos = pos + self.energy = energy + self.theta = float(theta) + self.phi = float(kwargs.pop('phi', 0.)) + self.depletion_thickness = float(kwargs.pop('depletion_thickness', 26.3)) + #todo: NOTE: in unit of um !!!!!!! + # Set other possible systematics with kwargs, including + # electric fields, pixel size, etc. + + self.path_to_root = os.getenv('DECO_ROOT_PATH') + self.path_to_geant = os.getenv('DECO_GEANT_PATH') + self.path_to_allpix = os.getenv('DECO_ALLPIX_PATH') + self.base_command = self.path_to_allpix + ' -c ./htc_wildfire/source_measurement.conf -o' + + def write_source_file(self, n_events, want_plot): + + with open('./htc_wildfire/source_measurement_replace.conf', 'r') as f: + data = f.readlines() + data[3] = data[3].format(n_events) + + data[18] = data[18].format(str(self.pos[0]) + " " + str(self.pos[1]) + " " + str(self.pos[2]) + "um") + + theta = pi - radians(self.theta) + phi = radians(self.phi) + dirx = 1*sin(theta)*cos(phi) + diry = 1*sin(theta)*sin(phi) + dirz = 1*cos(theta) + + data[20] = data[20].format(str(dirx) + " " + str(diry) + " " + str(dirz)) + + if want_plot == 'true': + data[38] = data[38].format('true' + "\n" + "output_linegraphs = true \n output_plots_step = 100ps \n output_plots_align_pixels = true \n output_plots_use_pixel_units = true") + else: + data[38] = data[38].format("false") + + + with open('./htc_wildfire/source_measurement.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + """ + def write_detector_file(self): + # USE THE PARAMETERS TO REWRITE DETECTOR FILE + with open('./htc_wildfire/detector_replace.conf', 'r') as f: + data = f.readlines() + data[-2] = data[-2].format(self.phi, self.theta, 0) + with open('./htc_wildfire/detector.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + with open('./htc_wildfire/htc_wildfire_shielded_replace.conf', 'r') as f: + data = f.readlines() + data[4] = data[4].format(self.depletion_thickness) + with open('./htc_wildfire/htc_wildfire_shielded.conf', 'w') as wf: + wf.writelines(data) + wf.close() + """ + + def set_output_file_name(self): + # write unique file name depending on parameters + if not os.path.exists("./output"): + os.system("mkdir ./output") + if not os.path.exists("./output/" + str(self.pid)): + os.system("mkdir ./output/" + str(self.pid)) + + self.outfile = str(self.pid) + "/" + "{}_theta_{:.1f}_phi_{:.1f}_thickiness_{:.1f}_highstats.txt".format(self.energy, self.theta, self.phi, self.depletion_thickness) + pass + + + def get_output_file_name(self): + try: + return self.outfile + except: + self.set_output_file_name() + return self.outfile + + + def run_simulation(self, n_events, want_charge_plot='false'): + + self.source_local_env() + # Run the allpix simulation + + output_file = self.get_output_file_name() + self.write_source_file(n_events, want_charge_plot) + + my_command = self.base_command[:] + + my_command += ' DepositionGeant4.particle_type="{}"'.format(self.pid) + my_command += ' -o DepositionGeant4.source_energy="{}"'.format(self.energy) + my_command += ' -o TextWriter.file_name="' + output_file + '"' + + """check if simulated using file name""" + if self.check_if_simulated(output_file) is True: + print("has been simulated") + return + + subprocess.call(my_command, shell=True) + + return + + + def check_if_simulated(self, filename): + # Check to see if simulation has already been run for + # this set of parameters + + if not os.path.exists("./output/" + filename): + return False + + return True + + + + def source_local_env(self): + geant = self.path_to_geant + "/bin/geant4.sh" + root = self.path_to_root + "/bin/thisroot.sh" + + pass + # Run source scripts for ROOT and GEANT if the + # simulation doesn't work + + +dep_thickness = 26.3 #um +pos = [0, 0, 0] + +energy = np.logspace(4, 5, 10) + +energy_list = [] + +for ene in energy: + energy_list.append(str(round(ene/pow(10, 6), 3))+"MeV") + + +#angles = ['0', '15', '30', '45', '60', '75'] + +#phis = ['0', '15', '30', '45', '60', '75', '90'] + +#energy = ['10GeV'] +angles = ['45'] +phis = ['0'] +particle_type = 'e-' + + +want_charge_plot = "false" + +for ene in energy_list: + for ang in angles: + for azi in phis: + a = DECOMuonSimulator(particle_type, ene, pos, ang, phi=azi, depletion_thickness=str(dep_thickness)) + + a.run_simulation(1000, want_charge_plot) + + diff --git a/Electron_energy_caliberator_analyze.py b/Electron_energy_caliberator_analyze.py new file mode 100644 index 0000000..b684ff2 --- /dev/null +++ b/Electron_energy_caliberator_analyze.py @@ -0,0 +1,206 @@ +import numpy as np +import pandas as pd +import matplotlib.pylab as plt +from matplotlib.colors import ListedColormap +import seaborn as sns +import os +plt.style.use("IceCube") +import matplotlib as mpl +from math import * +import sys + + + +def get_dEdx(M, E, Z, type): + coeff = 2.32 * 14. / 28.0855 * 0.307 + meev = .511 + gamma = E / M + 1 + beta = np.sqrt(1 - (1 / gamma ** 2)) + beta2gamma2 = (beta * gamma) ** 2 + if type != 'e+' and type != 'e-': + Tmax = (2 * meev * beta2gamma2) / (1 + 2 * gamma * meev / M + (meev / M) ** 2) + else: + Tmax = E + M + I = 10 * 14.0e-6 + hw = np.sqrt(2.32 * 14 / 28.0855) * 28.816e-6 + ln = np.log((2 * meev * beta2gamma2 * Tmax) / I ** 2) + delta_over_2 = np.log(hw / I) + np.log(np.sqrt(beta2gamma2)) - 0.5 + + return coeff * Z ** 2 / beta ** 2 * (0.5 * ln - beta ** 2 - delta_over_2) + + +def getdE(E_ini, dep_leng, p_type): + + E_curr = E_ini + remain_dis = dep_leng * pow(10, -6) + + while 1: + if E_curr <= 0 and remain_dis > 0: + return E_ini + + if remain_dis <= 0: + return E_ini - E_curr + + if E_curr < 1e5: + step_size = 0.001 * 1e-6 + elif E_curr > 1e5 and E_curr < 1e6: + step_size = 0.01 * 1e-6 + elif E_curr > 1e6 and E_curr < 5 * 1e6: + step_size = 0.1 * 1e-6 + elif E_curr > 5*1e6 and E_curr < 1e7: + step_size = 1e-6 + else: + step_size = 2 * 1e-6 + + dE_dx = get_dEdx(Mass, E_curr/pow(10, 6), 1.0, p_type) + + E_curr -= dE_dx*pow(10, 6) * step_size/0.01 + + remain_dis -= step_size + + + +class DECOLeptonAnalyzer(): + r'''This class is for making plots like the ones in the + notebook that read in a bunch of simulated files + and make analysis level plots''' + + def __init__(self, pid, phi, thickness): + self.pid = pid + + self.thichness = thickness + + self.phi = phi + + + def read_hit_file(self, filename): + f = open(filename, 'r') + + xhits, yhits, charge = [], [], [] + + # skip first 2 lines + f.readline() + f.readline() + + x, y, c = [], [], [] + + flag = 1 + while 1: + + temp = f.readline().split() + + if len(temp) < 1 or temp[0] == '#': + break + + if temp[0] == '===': + continue + if temp[0] == '---': + if flag == 1: + flag = 0 + else: + xhits.append(x) + yhits.append(y) + charge.append(c) + x, y, c = [], [], [] + + else: + x.append(float(temp[1][:-1])) + y.append(float(temp[2][:-1])) + c.append(float(temp[3][:-1])) + + if len(x) > 0: + xhits.append(x) + yhits.append(y) + charge.append(c) + + return xhits, yhits, charge + + + + + def get_E_deposited(self, en, ang): + x, y, c = self.read_hit_file( + "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(ang), + float(self.phi), self.thichness)) + + dE_list = [] + + for i in range(len(c)): + dE_list.append(np.array(c[i]).sum() * 3.62) + + + return np.array(dE_list) + + + + + +energy = np.logspace(4, 10, 100) + +energy_list = [] + +for ene in energy: + energy_list.append(str(ene/pow(10, 6))+"MeV") + +theta = '45' +phi = '0' +particle_type = 'mu+' + +thickness = 26.3 + +a = DECOLeptonAnalyzer(particle_type, phi, thickness) + +E_calibrate_list = [] +err_list = [] +E_exp_list = [] + +for i in range(len(energy)): + + print("working on " + str(round(energy[i]/pow(10, 6), 4)) + "MeV") + + dE_list = a.get_E_deposited(energy_list[i], theta) + + E_calibrate_list.append(dE_list.mean()) + + err_list.append(dE_list.std()) + + if particle_type == 'mu+' or particle_type == 'mu-': + Mass = 105.658 + if particle_type == 'e-' or particle_type == 'e+': + Mass = 0.511 + + E_exp_list.append(getdE(energy[i], thickness/2, particle_type)) + + + f = plt.figure() + + plt.hist(dE_list, bins=int((dE_list.max() - dE_list.min())/1000) + 1, label="initial E = " + str(round(energy[i]/pow(10, 6), 4)) + "MeV \nmean of calibrated E = " + + str(round(dE_list.mean()/pow(10, 6), 4)) + + "MeV \nstd = " + str(round(dE_list.std()/pow(10, 6), 4)) + "MeV") + + plt.xlabel("calibrated deposited energy in eV") + + plt.ylabel("# events (" + str(len(dE_list)) + " in total)") + + plt.legend() + + if not os.path.exists("./Individual_E_Calib_Plots"): + os.mkdir("./Individual_E_Calib_Plots") + if not os.path.exists("./Individual_E_Calib_Plots/" + str(particle_type)): + os.mkdir("./Individual_E_Calib_Plots/" + str(particle_type)) + + plt.savefig("./Individual_E_Calib_Plots/" + str(particle_type) + "/Eini_" + str(round(energy[i]/pow(10, 6), 4)) + "MeV.png", bbox_inches='tight') + + plt.close() + #plt.show() + + + +f = plt.figure() +plt.errorbar(energy, E_calibrate_list, yerr=err_list, label='Calibrated energy with uncertainty') +plt.plot(energy, E_exp_list, label='Expected fully reconstructed energy with uncertainty', linewidth=3) +plt.legend() +plt.xlabel("Initial energy in eV") +plt.ylabel("Calibrated energy in eV") +plt.xscale('log') +plt.show() diff --git a/LeptonAnalyzer.py b/LeptonAnalyzer.py index e055006..9a69d03 100644 --- a/LeptonAnalyzer.py +++ b/LeptonAnalyzer.py @@ -5,14 +5,14 @@ import seaborn as sns import os plt.style.use("IceCube") - +import matplotlib as mpl class DECOLeptonAnalyzer(): r'''This class is for making plots like the ones in the notebook that read in a bunch of simulated files and make analysis level plots''' - - def __init__(self, pid, energy_levels, en_float, theta_list, phi, thickness): + + def __init__(self, pid, energy_levels, en_float, theta_list, phi_list, thickness): self.pid = pid self.thichness = thickness @@ -24,9 +24,9 @@ def __init__(self, pid, energy_levels, en_float, theta_list, phi, thickness): self.en_float = en_float - self.theta_angles = theta_list + self.theta_list = theta_list - self.phi = phi + self.phi_list = phi_list self.data_list = self.data_processing() @@ -77,20 +77,23 @@ def read_hit_file(self, filename): return the distance by L2 norm of (x_max_diff, y_max_diff) """ def track_length(self, x, y): - x_dist_sq = np.power(np.max(x) - np.min(x), 2.) - y_dist_sq = np.power(np.max(y) - np.min(y), 2.) + x_dist_sq = np.power((np.max(x) - np.min(x)) * 0.9, 2.) + y_dist_sq = np.power((np.max(y) - np.min(y)) * 0.9, 2.) return np.power(x_dist_sq + y_dist_sq, 0.5) """ - + dE/dx in unit of MeV/cm """ - def get_dEdx(self, M, E, Z): + def get_dEdx(self, M, E, Z, type): coeff = 2.32 * 14. / 28.0855 * 0.307 meev = .511 gamma = E / M + 1 beta = np.sqrt(1 - (1 / gamma ** 2)) beta2gamma2 = (beta * gamma) ** 2 - Tmax = (2 * meev * beta2gamma2) / (1 + 2 * gamma * meev / M + (meev / M) ** 2) + if type != 'e+' and type != 'e-': + Tmax = (2 * meev * beta2gamma2) / (1 + 2 * gamma * meev / M + (meev / M) ** 2) + else: + Tmax = E + M I = 10 * 14.0e-6 hw = np.sqrt(2.32 * 14 / 28.0855) * 28.816e-6 ln = np.log((2 * meev * beta2gamma2 * Tmax) / I ** 2) @@ -109,22 +112,88 @@ def data_processing(self): counter = 0 curr_energy = 0 for en in self.energy_levels: - for ang in self.theta_angles: - x, y, c = self.read_hit_file("./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(ang), float(self.phi), self.thichness)) - try: - for j in range(len(x)): - charge = np.sum(c[j]) - length = self.track_length(x[j], y[j]) - dE_dX = charge / np.power(length ** 2 + (self.thichness / 0.9) ** 2, 0.5) - data_list.loc[counter] = [en, ang, '30', charge, self.en_float[curr_energy], length, dE_dX] - counter += 1 - except: - continue + for theta in self.theta_list: + for phi in self.phi_list: + x, y, c = self.read_hit_file("./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(theta), float(phi), self.thichness)) + try: + for j in range(len(x)): + charge = np.sum(c[j]) + length = self.track_length(x[j], y[j]) + dQ_dX = charge / np.power(length ** 2 + self.thichness ** 2, 0.5) + data_list.loc[counter] = [en, theta, phi, charge, self.en_float[curr_energy], length, dQ_dX] + counter += 1 + except: + continue curr_energy += 1 return data_list + + def plot_single(self, en, theta, phi): + x, y, c = self.read_hit_file( + "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(theta), + float(phi), self.thichness)) + + for j in range(50): + # try: + image = np.zeros((4500, 4500)) + for i in range(len(x[j])): + image[int(y[j][i]), int(x[j][i])] = c[j][i] + + + med_x = np.median(x[j]) + med_y = np.median(y[j]) + size = 50. + + title = "1 GeV $\mu^{+}$" + fig1 = plt.figure(1, figsize=(8, 8)) + ax = fig1.add_subplot(111) + fig1.set_facecolor('white') + + # my_cmap = ListedColormap(sns.color_palette("Blues", 50)) + # my_cmap = ListedColormap(sns.palplot(sns.cubehelix_palette(8, start=2, rot=0, dark=0, light=.95, reverse=True))) + my_cmap = mpl.cm.hot + + # image = np.where(image == 0.0, np.nan, image) + + im = ax.imshow(image, cmap=my_cmap, # interpolation="gaussian", + aspect="auto", vmax=100., vmin=0.0) + ax.set_xlim([med_x - size, med_x + size]) + ax.set_ylim([med_y - size, med_y + size]) + ax.set_xlabel("X (pixels)") + ax.set_ylabel("Y (pixels)") + ax.set_title(title) + + ax.grid(color="#ffffff") + cb = fig1.colorbar(im, orientation="vertical", + shrink=0.8, + fraction=0.05, + pad=0.15) + label = "Pixel Luminance" + cb.set_label(label) + ax.text(med_x + size * 0.3, med_y + size * 0.7, + "Simulation", fontsize=24, color='w', weight='heavy') + plt.show() + def bethe_bloch_plot(self): + me = .511 + E_array = np.logspace(-2., 4., 100) + BB = [] + + for E in E_array: + + if self.pid == 'mu+' or self.pid == 'mu-': + BB.append(self.get_dEdx(206.7 * me, E, 1., self.pid)) + elif self.pid == 'e+' or self.pid == 'e-': + BB.append(self.get_dEdx(me, E, 1., self.pid)) + + + + BB = np.array(BB) + + #todo: may change this /2 if counting hole charge is possible + BB = BB * 1e6 / 2 * (1. / 3.62) * 1e-4 # convert from MeV / cm to electron charge per um + fig, ax = plt.subplots(figsize=(9, 6)) @@ -135,12 +204,12 @@ def bethe_bloch_plot(self): h = plt.hist2d(np.log10(self.data_list['Energy (GeV)']), np.log10(self.data_list['Charge per unit length']), bins=[np.linspace(-1., 5., 14), np.linspace(1.5, 3.5, 35)], cmin=1., cmap=my_cmap) plt.colorbar(label="Number of Events") - plt.title("Muon Losses") - #plt.plot(np.log10(E_array), np.log10(muon_BB), c='r', label="Bethe-Bloch", lw=3) + plt.title(str(self.pid) + " Losses") + plt.plot(np.log10(E_array), np.log10(BB), c='r', label="Bethe-Bloch", lw=3) plt.ylabel(r'$\log (\frac{dQ}{dx} \times \frac{0.9 \mu m}{q_{e}})$', fontsize=26) - plt.xlabel('$\log$( $E_{\mu}$ / MeV) ') - plt.xlim(0, 4.2) - plt.ylim(1.5, 3.25) + plt.xlabel('$\log$( $E_{' + str(self.pid) + '}$ / MeV) ') + #plt.xlim(-2, 4.2) + #plt.ylim(1.5, 3.25) plt.show() @@ -148,13 +217,14 @@ def check_if_all_simulated(self): missing_list = [] for en in self.energy_levels: - for ang in self.theta_angles: + for theta in self.theta_list: + for phi in self.phi_list: - curr_file_name = "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format( - self.pid, en, float(ang), float(phi), self.thichness) + curr_file_name = "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format( + self.pid, en, float(theta), float(phi), self.thichness) - if not os.path.exists(curr_file_name): - missing_list.append(curr_file_name) + if not os.path.exists(curr_file_name): + missing_list.append(curr_file_name) if len(missing_list) == 0: print("find all required files:") @@ -182,7 +252,7 @@ def track_length_vs_angle_violinplot(self): index = 0 for i in range(self.data_list.__len__()): if self.data_list.loc[i]['Energy'] == '100MeV' or self.data_list.loc[i]['Energy'] == '10GeV': - if self.data_list.loc[i]['Track Length (pixels)'] <= 250 and self.data_list.loc[i]['Track Length (pixels)'] > 0: + if self.data_list.loc[i]['Track Length (pixels)'] <= 250: curr_data.loc[index] = self.data_list.loc[i] index += 1 @@ -215,13 +285,17 @@ def hillas_width_histogram(self): theta_list = ['0', '15', '30', '45', '60', '75'] -phi = 30 +#phi_list = ['0', '15', '30', '45', '60', '75', '90'] + +phi_list = ['0'] thickness = 26.3 -a = DECOLeptonAnalyzer('mu+', energy_levels, en_float, theta_list, phi, thickness) +a = DECOLeptonAnalyzer('e-', energy_levels, en_float, theta_list, phi_list, thickness) + +#a.track_length_vs_angle_violinplot() +a.bethe_bloch_plot() -a.track_length_vs_angle_violinplot() -#a.bethe_bloch_plot() +#a.plot_single('1GeV', '75') #x, y, c = a.read_hit_file('./output/mu+/100KeV_theta_45.0_phi_30.0_thickiness_26.3_highstats.txt') diff --git a/MuonSimulator.py b/MuonSimulator.py index ddc4b23..6154ae8 100644 --- a/MuonSimulator.py +++ b/MuonSimulator.py @@ -3,13 +3,15 @@ import os, sys import subprocess from pipes import quote +from math import * class DECOMuonSimulator(): r'''Simulator class for muons that uses allpix and GEANT''' - def __init__(self, pid, energy, theta, **kwargs): + def __init__(self, pid, energy, pos, theta, **kwargs): self.pid = pid #Particle id, ie mu+, e- + self.pos = pos self.energy = energy self.theta = float(theta) self.phi = float(kwargs.pop('phi', 0.)) @@ -23,16 +25,33 @@ def __init__(self, pid, energy, theta, **kwargs): self.path_to_allpix = os.getenv('DECO_ALLPIX_PATH') self.base_command = self.path_to_allpix + ' -c ./htc_wildfire/source_measurement.conf -o' - def write_source_file(self, n_events): + def write_source_file(self, n_events, want_plot): with open('./htc_wildfire/source_measurement_replace.conf', 'r') as f: data = f.readlines() data[3] = data[3].format(n_events) + + data[18] = data[18].format(str(self.pos[0]) + " " + str(self.pos[1]) + " " + str(self.pos[2]) + "um") + + theta = pi - radians(self.theta) + phi = radians(self.phi) + dirx = 1*sin(theta)*cos(phi) + diry = 1*sin(theta)*sin(phi) + dirz = 1*cos(theta) + + data[20] = data[20].format(str(dirx) + " " + str(diry) + " " + str(dirz)) + + if want_plot == 'true': + data[38] = data[38].format('true' + "\n" + "output_linegraphs = true \n output_plots_step = 100ps \n output_plots_align_pixels = true \n output_plots_use_pixel_units = true") + else: + data[38] = data[38].format("false") + + with open('./htc_wildfire/source_measurement.conf', 'w') as wf: wf.writelines(data) wf.close() - + """ def write_detector_file(self): # USE THE PARAMETERS TO REWRITE DETECTOR FILE with open('./htc_wildfire/detector_replace.conf', 'r') as f: @@ -48,7 +67,7 @@ def write_detector_file(self): with open('./htc_wildfire/htc_wildfire_shielded.conf', 'w') as wf: wf.writelines(data) wf.close() - + """ def set_output_file_name(self): # write unique file name depending on parameters @@ -69,14 +88,13 @@ def get_output_file_name(self): return self.outfile - def run_simulation(self, n_events=100): + def run_simulation(self, n_events, want_charge_plot='false'): self.source_local_env() # Run the allpix simulation output_file = self.get_output_file_name() - self.write_detector_file() - self.write_source_file(n_events) + self.write_source_file(n_events, want_charge_plot) my_command = self.base_command[:] @@ -114,17 +132,28 @@ def source_local_env(self): # simulation doesn't work +dep_thickness = 26.3 #um +pos = [0, 0, dep_thickness/2] - -energy = ['10keV', '31.6keV', '100keV', '316keV', '1MeV', '3.16MeV', - '10MeV', '31.6MeV', '100MeV', '316MeV', '1GeV', '3.16GeV', '10GeV'] +energy = ['10keV', '31.6keV', '100keV', '316keV', '1MeV', '3.16MeV', '10MeV', '31.6MeV', '100MeV', '316MeV', '1GeV', '3.16GeV', '10GeV'] angles = ['0', '15', '30', '45', '60', '75'] -for e in energy: +#phis = ['0', '15', '30', '45', '60', '75', '90'] + +#energy = ['10GeV'] +#angles = ['45'] +phis = ['0'] +particle_type = 'e-' + + +want_charge_plot = "false" + +for ene in energy: for ang in angles: - a = DECOMuonSimulator('mu+', e, ang, phi='30', depletion_thickness='26.3') + for azi in phis: + a = DECOMuonSimulator(particle_type, ene, pos, ang, phi=azi, depletion_thickness=str(dep_thickness)) - a.run_simulation(100) + a.run_simulation(100, want_charge_plot) diff --git a/Photon_corss_section_Simulator.py b/Photon_corss_section_Simulator.py new file mode 100644 index 0000000..40fe778 --- /dev/null +++ b/Photon_corss_section_Simulator.py @@ -0,0 +1,168 @@ +r'''Class that handles various inputs and systematic parameters + and simulates muons''' +import os, sys +import subprocess +from pipes import quote +from math import * +import numpy as np + +class DECOMuonSimulator(): + r'''Simulator class for muons that uses allpix + and GEANT''' + + def __init__(self, pid, energy, pos, theta, **kwargs): + self.pid = pid #Particle id, ie mu+, e- + self.pos = pos + self.energy = energy + self.theta = float(theta) + self.phi = float(kwargs.pop('phi', 0.)) + self.depletion_thickness = float(kwargs.pop('depletion_thickness', 26.3)) + #todo: NOTE: in unit of um !!!!!!! + # Set other possible systematics with kwargs, including + # electric fields, pixel size, etc. + + self.path_to_root = os.getenv('DECO_ROOT_PATH') + self.path_to_geant = os.getenv('DECO_GEANT_PATH') + self.path_to_allpix = os.getenv('DECO_ALLPIX_PATH') + self.base_command = self.path_to_allpix + ' -c ./htc_wildfire/source_measurement.conf -o' + + def write_source_file(self, n_events, want_plot): + + with open('./htc_wildfire/source_measurement_replace.conf', 'r') as f: + data = f.readlines() + data[3] = data[3].format(n_events) + + data[18] = data[18].format(str(self.pos[0]) + " " + str(self.pos[1]) + " " + str(self.pos[2]) + "um") + + theta = pi - radians(self.theta) + phi = radians(self.phi) + dirx = 1*sin(theta)*cos(phi) + diry = 1*sin(theta)*sin(phi) + dirz = 1*cos(theta) + + data[20] = data[20].format(str(dirx) + " " + str(diry) + " " + str(dirz)) + + if want_plot == 'true': + data[38] = data[38].format('true' + "\n" + "output_linegraphs = true \n output_plots_step = 100ps \n output_plots_align_pixels = true \n output_plots_use_pixel_units = true") + else: + data[38] = data[38].format("false") + + + with open('./htc_wildfire/source_measurement.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + """ + def write_detector_file(self): + # USE THE PARAMETERS TO REWRITE DETECTOR FILE + with open('./htc_wildfire/detector_replace.conf', 'r') as f: + data = f.readlines() + data[-2] = data[-2].format(self.phi, self.theta, 0) + with open('./htc_wildfire/detector.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + with open('./htc_wildfire/htc_wildfire_shielded_replace.conf', 'r') as f: + data = f.readlines() + data[4] = data[4].format(self.depletion_thickness) + with open('./htc_wildfire/htc_wildfire_shielded.conf', 'w') as wf: + wf.writelines(data) + wf.close() + """ + + def set_output_file_name(self): + # write unique file name depending on parameters + if not os.path.exists("./output"): + os.system("mkdir ./output") + if not os.path.exists("./output/" + str(self.pid)): + os.system("mkdir ./output/" + str(self.pid)) + + self.outfile = str(self.pid) + "/" + "{}_theta_{:.1f}_phi_{:.1f}_thickiness_{:.1f}_highstats.txt".format(self.energy, self.theta, self.phi, self.depletion_thickness) + pass + + + def get_output_file_name(self): + try: + return self.outfile + except: + self.set_output_file_name() + return self.outfile + + + def run_simulation(self, n_events, want_charge_plot='false'): + + self.source_local_env() + # Run the allpix simulation + + output_file = self.get_output_file_name() + self.write_source_file(n_events, want_charge_plot) + + my_command = self.base_command[:] + + my_command += ' DepositionGeant4.particle_type="{}"'.format(self.pid) + my_command += ' -o DepositionGeant4.source_energy="{}"'.format(self.energy) + my_command += ' -o TextWriter.file_name="' + output_file + '"' + + """check if simulated using file name""" + if self.check_if_simulated(output_file) is True: + print("has been simulated") + return + + subprocess.call(my_command, shell=True) + + return + + + def check_if_simulated(self, filename): + # Check to see if simulation has already been run for + # this set of parameters + + if not os.path.exists("./output/" + filename): + return False + + return True + + + + def source_local_env(self): + geant = self.path_to_geant + "/bin/geant4.sh" + root = self.path_to_root + "/bin/thisroot.sh" + + pass + # Run source scripts for ROOT and GEANT if the + # simulation doesn't work + + +dep_thickness = 26.3 #um +pos = [0, 0, dep_thickness/2] + +#energy = ['10keV', '31.6keV', '100keV', '316keV', '1MeV', '3.16MeV', '10MeV', '31.6MeV', '100MeV', '316MeV', '1GeV', '3.16GeV', '10GeV'] + +#angles = ['0', '15', '30', '45', '60', '75'] + +#phis = ['0', '15', '30', '45', '60', '75', '90'] + +#energy = ['10GeV'] + +energy = np.logspace(4, 9, 100) + +energy_list = [] + +for ene in energy: + energy_list.append(str(round(ene/pow(10, 6), 3))+"MeV") + +angles = ['45'] +phis = ['0'] +particle_type = 'gamma' + + +want_charge_plot = "false" + +for ene in energy_list: + for ang in angles: + for azi in phis: + a = DECOMuonSimulator(particle_type, ene, pos, ang, phi=azi, depletion_thickness=str(dep_thickness)) + + a.run_simulation(20000, want_charge_plot) + + diff --git a/Photon_cross_section_analyer.py b/Photon_cross_section_analyer.py new file mode 100644 index 0000000..9465c8e --- /dev/null +++ b/Photon_cross_section_analyer.py @@ -0,0 +1,117 @@ +import numpy as np +import pandas as pd +import matplotlib.pylab as plt +from matplotlib.colors import ListedColormap +import seaborn as sns +import os +plt.style.use("IceCube") +import matplotlib as mpl + +class DECOLeptonAnalyzer(): + r'''This class is for making plots like the ones in the + notebook that read in a bunch of simulated files + and make analysis level plots''' + + def __init__(self, pid, phi, thickness): + self.pid = pid + + self.thichness = thickness + + self.phi = phi + + + def read_hit_file(self, filename): + f = open(filename, 'r') + + xhits, yhits, charge = [], [], [] + + # skip first 2 lines + f.readline() + f.readline() + + x, y, c = [], [], [] + + flag = 1 + while 1: + + temp = f.readline().split() + + if len(temp) < 1 or temp[0] == '#': + break + + if temp[0] == '===': + continue + if temp[0] == '---': + if flag == 1: + flag = 0 + else: + xhits.append(x) + yhits.append(y) + charge.append(c) + x, y, c = [], [], [] + + else: + x.append(float(temp[1][:-1])) + y.append(float(temp[2][:-1])) + c.append(float(temp[3][:-1])) + + if len(x) > 0: + xhits.append(x) + yhits.append(y) + charge.append(c) + + return xhits, yhits, charge + + + + + def get_events_num(self, en, ang): + x, y, c = self.read_hit_file( + "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(ang), + float(self.phi), self.thichness)) + + return len(x) + + def plot_cross_section(self, ene_number_list, ene_list, theta): + + num_events = [] + + for ene in ene_list: + + num = self.get_events_num(ene, theta) + + num_events.append(num) + + num_event_err = np.sqrt(np.array(num_events)) + + f = plt.figure() + + plt.errorbar(ene_number_list, num_events, yerr=num_event_err) + + plt.xscale('log') + plt.yscale('log') + + plt.xlabel('E in eV') + plt.ylabel('log(#Events)') + + plt.show() + + + +energy = np.logspace(4, 9, 100) + +energy_list = [] + +for ene in energy: + energy_list.append(str(round(ene/pow(10, 6), 3))+"MeV") + +theta = '45' +phi = '0' +particle_type = 'gamma' + +thickness = 26.3 + +a = DECOLeptonAnalyzer(particle_type, phi, thickness) + +a.plot_cross_section(energy, energy_list, theta) + diff --git a/README.md b/README.md index 3894402..0f53323 100644 --- a/README.md +++ b/README.md @@ -1 +1,3 @@ -This repository is for the simulation of cell phone camera image sensors and relies on the installation of both Allpix-Squared and GEANT4. +This repository is for the simulation of cell phone camera image sensors and relies on the installation of both Allpix-Squared and GEANT4. + +In Simulation_with_CNN folder, a CNN is used to do event classification after simulation after some post image processings. diff --git a/Simultaion_with_CNN/CNN_on_simulation/CNN_test.py b/Simultaion_with_CNN/CNN_on_simulation/CNN_test.py new file mode 100644 index 0000000..87b6a84 --- /dev/null +++ b/Simultaion_with_CNN/CNN_on_simulation/CNN_test.py @@ -0,0 +1,133 @@ +import json +from PIL import Image +import argparse +import os +import h5py +import zoom_predictions_4class +import numpy as np +import pandas as pd +import decotools as dt +from zoom_predictions_4class import run_blob_classifier +from shutil import copy2 + + + +# CNN class probability thresholds +# Event is classified as 'key' if CNN prob is > value, otherwise 'ambiguous' +class_thresholds = {'track': 0.8, + 'worm': 0.8, + 'spot': 0.8, + 'noise': 0.8 + } + +def get_predicted_label(event, thresholds=class_thresholds): + """Return predicted label for an event. + + Parameters + ---------- + event : pandas.Series + Row of dataframe containing event id, image path, centroid coordinates, + and CNN class probabilities + thresholds : dict + Dictionary containing cutoff threshold for each class probability + e.g., 'track' : 0.8 + + Returns + ------- + pred : str + Predicted event classification label + """ + + if event['p_track'] < 0: + pred = 'Edge' + elif event['p_noise'] > thresholds['noise']: + pred = 'Noise' + elif event['p_spot'] > thresholds['spot']: + pred = 'Spot' + elif event['p_track'] > thresholds['track']: + pred = 'Track' + elif event['p_worm'] > thresholds['worm']: + pred = 'Worm' + else: + pred = 'Ambiguous' + + return pred + + +def get_events(): + + # Replace extension with the jpg extension of the raw, i.e. unprocessed, event image + absolute_paths = [] + + + base = '../MC_Simulation/output/' + file_dir + '/' + source_dir + '/' + + for filename in os.listdir(base): + if filename.endswith(".txt"): + absolute_paths.append(base + filename) + + return absolute_paths + +file_dir = "qdc_res_0_smear_0_slope_7" + +source_dir = 'mu+_2Dimages_fac1.5_dev0_condev0' + +if __name__ == '__main__': + + cnn_weights_file = './final_trained_weights.h5' + + # Get all paths + paths = get_events() + print(paths) + # Use args.index and args.n_images for batch processing + df = run_blob_classifier(paths, 4, + weights_file=cnn_weights_file) + + if not os.path.exists("./classification"): + os.mkdir("./classification") + + if not os.path.exists("./classification/" + file_dir): + os.mkdir("./classification/" + file_dir) + + if not os.path.exists("./classification/" + file_dir + '/' + source_dir): + os.mkdir("./classification/" + file_dir + '/' + source_dir) + + total_num = 0 + track_num = 0 + + for i in range(len(df["p_track"])): + + total_num += 1 + + if df["p_track"][i] > 0.5: + track_num += 1 + + if not os.path.exists("./classification/" + file_dir + '/' + source_dir + "/track"): + os.mkdir("./classification/" + file_dir + '/' + source_dir + "/track") + + copy2(paths[i], "./classification/" + file_dir + '/' + source_dir + "/track/") + + if df["p_worm"][i] > 0.5: + + if not os.path.exists("./classification/" + file_dir + '/' + source_dir + "/worm"): + os.mkdir("./classification/" + file_dir + '/' + source_dir + "/worm") + + copy2(paths[i], "./classification/" + file_dir + '/' + source_dir + "/worm/") + + if df["p_spot"][i] > 0.5: + + if not os.path.exists("./classification/" + file_dir + '/' + source_dir + "/spot"): + os.mkdir("./classification/" + file_dir + '/' + source_dir + "/spot") + + copy2(paths[i], "./classification/" + file_dir + '/' + source_dir + "/spot/") + + + print("track rate " + str(float(track_num)/total_num)) + + + + + + + + diff --git a/Simultaion_with_CNN/CNN_on_simulation/readme.txt b/Simultaion_with_CNN/CNN_on_simulation/readme.txt new file mode 100644 index 0000000..095c31e --- /dev/null +++ b/Simultaion_with_CNN/CNN_on_simulation/readme.txt @@ -0,0 +1,2 @@ +Apply the CNN used in DECO on simulated data to select tracks with P > 0.5 +The weight file can be found in https://www.dropbox.com/s/av58b92f0mvxtp6/final_trained_weights.h5?dl=0 diff --git a/Simultaion_with_CNN/CNN_on_simulation/zoom_predictions_4class.py b/Simultaion_with_CNN/CNN_on_simulation/zoom_predictions_4class.py new file mode 100644 index 0000000..82bae48 --- /dev/null +++ b/Simultaion_with_CNN/CNN_on_simulation/zoom_predictions_4class.py @@ -0,0 +1,185 @@ + +#------------------------------------------------------------------------# +# Author: Miles Winter # +# Date: 07-14-2017 # +# Project: DECO # +# Desc: zoom in on brightest pixel and classify blobs with CNN # +# Note: Need the following installed: # +# $ pip install --user --upgrade h5py theano keras # +# Change keras backend to theano (default is tensorflow) # +# Importing keras generates a .json config file # +# $ KERAS_BACKEND=theano python -c "from keras import backend" # +# Next, to change "backend": "tensorflow" -> "theano" type # +# $ sed -i 's/tensorflow/theano/g' $HOME/.keras/keras.json # +# Documentation at https://keras.io/backend/ # +#------------------------------------------------------------------------# + + +import os +import numpy as np +import pandas as pd +import keras +from keras.models import load_model +from keras.models import Sequential +from keras.layers import Conv2D, MaxPooling2D, Cropping2D +from keras.layers import Flatten, Dense, Dropout +from keras.layers.advanced_activations import LeakyReLU +from PIL import Image +from collections import defaultdict + + +def build_model(n_classes, training=False): + """Define model structure""" + model = Sequential() + + if training: + model.add(Cropping2D(cropping=18,input_shape=(64, 64, 1))) + model.add(Conv2D(64, (3, 3), padding='same')) + else: + model.add(Conv2D(64, (3, 3), padding='same', input_shape=(64, 64, 1))) + + model.add(LeakyReLU(alpha=0.3)) + model.add(Conv2D(64, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(MaxPooling2D(pool_size=(2, 2))) + model.add(Dropout(0.2)) + + model.add(Conv2D(128, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(Conv2D(128, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(MaxPooling2D(pool_size=(2, 2))) + model.add(Dropout(0.2)) + + model.add(Conv2D(256, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(Conv2D(256, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(MaxPooling2D(pool_size=(2, 2))) + model.add(Dropout(0.2)) + + model.add(Conv2D(512, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(Conv2D(512, (3, 3), padding='same')) + model.add(LeakyReLU(alpha=0.3)) + model.add(MaxPooling2D(pool_size=(2, 2))) + model.add(Dropout(0.2)) + + model.add(Flatten()) + model.add(Dense(2048)) + model.add(LeakyReLU(alpha=0.3)) + model.add(Dropout(0.5)) + model.add(Dense(2048)) + model.add(LeakyReLU(alpha=0.3)) + model.add(Dropout(0.5)) + model.add(Dense(n_classes, activation='softmax')) + return model + +def get_event_id(path): + """Return the event ID for a given image file.""" + directory, filename = os.path.split(path) + event_id, ext = os.path.splitext(filename) + return event_id + + +def get_predicted_label(probs, track_thresh): + """Returns predicted label. Track if prediction is > track_thresh """ + if probs[-1] >= track_thresh: + return 'Track' + else: + return 'Other' + + +def get_crop_range(maxX, maxY, size=32): + """define region of image to crop""" + return maxX - size, maxX + size, maxY - size, maxY + size + + +def pass_edge_check(maxX, maxY, img_shape, crop_size=64): + """checks if image is on the edge of the sensor""" + x0, x1, y0, y1 = get_crop_range(maxX, maxY, size=crop_size / 2) + checks = np.array([x0 >= 0, x1 <= img_shape[0], y0 >= 0, y1 <= img_shape[1]]) + return checks.all() == True + + +def convert_image(img, dims=64, channels=1): + """convert image to grayscale, normalize, and reshape""" + img = np.array(img, dtype='float32') + gray_norm_img = np.mean(img, axis=-1) + return gray_norm_img.reshape(1, dims, dims, channels) + + +def get_brightest_pixel(img): + """get brightest image pixel indices""" + img = np.array(img) + #print(img.max()) + return np.unravel_index(img.argmax(), img.shape) + + +def run_blob_classifier(paths, n_classes, track_thresh=0.8, + weights_file='/data/user/mrmeehan/deco/classification/cnn/final_trained_weights.h5'): + """Classify blobs with CNN""" + # Build CNN model and load weights + try: + model = build_model(n_classes) + model.load_weights(weights_file) + except IOError: + print('Weights could not be found ... check path') + raise SystemExit + + class_labels = ['worm', 'spot', 'track', 'noise'] + data = defaultdict(list) + + file_counter = 0 + + for filename in paths: + + print("reading " + str(file_counter) + " out of " + str(len(paths))) + file_counter += 1 + + f_reader = open(filename, 'r') + image = [] + while 1: + temp = f_reader.readline().split() + temp_hold = [] + if len(temp) < 1: + break + + for content in temp: + temp_hold.append(float(content)) + + image.append(temp_hold) + + image = np.array(image) + + # Find the brightest pixel + + predicted_label = '' + probability = 0. + # Check if blob is near sensor edge + + + # Convert to grayscale, normalize, and reshape + + gray_image = (np.array(image, dtype='float32')/255).reshape(1, 64, 64, 1) + + + # Predict image classification + probability = model.predict(gray_image, batch_size=1, verbose=0) + probability = probability.reshape(n_classes,) + + # Convert prediction probability to a single label + #predicted_label = get_predicted_label(probability, track_thresh) + + + # Add predicted class label and probabilities from model + for idx, class_label in enumerate(class_labels): + data['p_{}'.format(class_label)].append(probability[idx]) + + data['image_file'].append(filename) + event_id = get_event_id(filename) + data['event_id'].append(event_id) + + df_data = pd.DataFrame.from_records(data) + + return df_data diff --git a/Simultaion_with_CNN/MC_simulation/DECO_MC_simulator.py b/Simultaion_with_CNN/MC_simulation/DECO_MC_simulator.py new file mode 100644 index 0000000..8917f56 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/DECO_MC_simulator.py @@ -0,0 +1,38 @@ +import numpy as np +from math import * +from Simulator_interface import DECOSimulator +from random_muon_generator import get_rand_ini_muon +from random_e_gamma_generator import get_rand_ini_e_gamma + +# energy range and particle type setup +E_min = 0.01 # GeV +E_max = 100 # GeV +dep_thickness = 26.3 #um +pos = [0, 0, dep_thickness/2] +particle_type = 'mu+' + +# total number of events to be simulated +num = 2000 + +for i in range(num): + print("working on " + str(i) + " out of " + str(num)) + + # get a random muon. Can also use electron or muon generator to generate isotropic zenith angle and uniform energy + ene, zen, azi = get_rand_ini_muon(E_min, E_max) + + zenith = str(degrees(zen)) + azimuth = str(degrees(azi)) + + energy = str(ene) + "GeV" + + want_charge_plot = "false" + + # qdc parameters, electronic noise, bias voltage, gain factor, pixel size, and seed (random if not specified) + a = DECOSimulator(particle_type, energy, pos, zenith, phi=azimuth, + depletion_thickness=str(dep_thickness), qdc_resolution=0, qdc_smearing=0, qdc_slope=0, + electronics_noise=0, bias_v=-25, gain=1.0, pixel_size=0.9, seed=10086) + + # run one event fer this setup + a.run_simulation(1, want_charge_plot) + + diff --git a/Simultaion_with_CNN/MC_simulation/Image_process_colored.py b/Simultaion_with_CNN/MC_simulation/Image_process_colored.py new file mode 100644 index 0000000..c1ce199 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/Image_process_colored.py @@ -0,0 +1,321 @@ +import numpy as np +import os +from math import * +from skimage import * +from PIL import Image, ImageEnhance + +""" +for HTC A510 we have a 2592x1944 CMOS sensor, we assign Bayer filter as: +green red +blue green +so for a pixel at (x, y), if x % 2 == 0 and y % 2 == 0 it is blue, if x % 2 == 1 and y % 2 == 1 it is red +""" +def getlabel(posx, posy): + + if posx % 2 == 0 and posy % 2 == 0: + return 'b' + elif posx % 2 == 1 and posy % 2 == 1: + return 'r' + else: + return 'g' + +"""get the background data""" +def getdata(filename): + f = open(filename, 'r') + + output = [] + + while 1: + temp = f.readline().split() + + output_temp = [] + + if len(temp) < 1: + break + + for num in temp: + output_temp.append(float(num)) + + output.append(output_temp) + + return np.array(output) + + +"""check if the event is on the edge of image""" +def pass_edge_check(maxX, maxY, img_shape, crop_size=64): + """checks if image is on the edge of the sensor""" + x0, x1, y0, y1 = get_crop_range(maxX, maxY, size=crop_size / 2) + checks = np.array([x0 >= 0, x1 <= img_shape[0], y0 >= 0, y1 <= img_shape[1]]) + return checks.all() == True + + +"""check the range to crop 64x64 image""" +def get_crop_range(maxX, maxY, size=32): + """define region of image to crop""" + return maxX - size, maxX + size, maxY - size, maxY + size + + +"""find the brightest pixel, since image crop is centered at the brightest point""" +def get_brightest_pixel(img): + """get brightest image pixel indices""" + img = np.array(img) + return np.unravel_index(img.argmax(), img.shape) + + +"""read the output file from simulation""" +def read_hit_file(filename): + f = open(filename, 'r') + + xhits, yhits, charge = [], [], [] + + # skip first 2 lines + f.readline() + f.readline() + + x, y, c = [], [], [] + + flag = 1 + while 1: + + temp = f.readline().split() + + if len(temp) < 1 or temp[0] == '#': + break + + if temp[0] == '===': + continue + if temp[0] == '---': + if flag == 1: + flag = 0 + else: + xhits.append(x) + yhits.append(y) + charge.append(c) + x, y, c = [], [], [] + + else: + x.append(float(temp[1][:-1])) + y.append(float(temp[2][:-1])) + c.append(float(temp[3][:-1])) + + if len(x) > 0: + xhits.append(x) + yhits.append(y) + charge.append(c) + + return xhits, yhits, charge + + +"""Color interpolation using mean algorithm and Bayer filter""" +def color_processing(image_arr): + + new_img = np.zeros((64, 64, 3)) + + for i in range(len(image_arr)): + for j in range(len(image_arr[i])): + + """do not work on edge pixels""" + if i == 0 or i == len(image_arr) - 1: + continue + if j == 0 or j == len(image_arr[i]) - 1: + continue + + label = getlabel(j, i) + + if label == 'r': + r = image_arr[i][j] + g = float(image_arr[i - 1][j] + image_arr[i + 1][j] + image_arr[i][j - 1] + image_arr[i][j + 1])/4 + b = float(image_arr[i-1][j-1] + image_arr[i + 1][j + 1] + image_arr[i + 1][j - 1] + image_arr[i - 1][j + 1])/4 + + elif label == 'b': + r = float(image_arr[i-1][j-1] + image_arr[i + 1][j + 1] + image_arr[i + 1][j - 1] + image_arr[i - 1][j + 1])/4 + g = float(image_arr[i - 1][j] + image_arr[i + 1][j] + image_arr[i][j - 1] + image_arr[i][j + 1])/4 + b = image_arr[i][j] + + + elif label == 'g' and getlabel(j - 1, i) == 'r': + + r = float(image_arr[i][j - 1] + image_arr[i][j + 1]) / 2 + + g = image_arr[i][j] + + b = float(image_arr[i - 1][j] + image_arr[i + 1][j]) / 2 + + + elif label == 'g' and getlabel(j - 1, i) == 'b': + + b = float(image_arr[i][j - 1] + image_arr[i][j + 1]) / 2 + + g = image_arr[i][j] + + r = float(image_arr[i - 1][j] + image_arr[i + 1][j]) / 2 + + else: + r, g, b = -1, -1, -1 + print("error!!!") + exit() + + # if r/255 < 0.0031308: + # r = r * 12.92 + # else: + # r = (1.055 * pow(r/255, 1.0/2.4) - 0.055)*255 + # + # + # if g/255 < 0.0031308: + # g = g * 12.92 + # else: + # g = (1.055 * pow(g/255, 1.0/2.4) - 0.055)*255 + # + # + # if b/255 < 0.0031308: + # b = b * 12.92 + # else: + # b = (1.055 * pow(b/255, 1.0/2.4) - 0.055)*255 + + new_img[i][j][0] = r #0.299 * r + 0.587 * g + 0.114 * b + new_img[i][j][1] = g + new_img[i][j][2] = b + + return new_img + + +"""Add background from data events""" +def background_addition(image_arr): + + + back_num = len([f for f in os.listdir("../CNN_on_training_data/back_by_device/" + model + "/") if f.endswith('.txt')]) + + index_back = np.random.randint(0, back_num - 1) + + back_arr = getdata("../CNN_on_training_data/back_by_device/" + model + "/" + str(index_back) + ".txt") + + image_arr += back_arr + + for i in range(len(image_arr)): + for j in range(len(image_arr[i])): + if image_arr[i][j] > 255: + image_arr[i][j] = 255 + image_arr[i][j] = int(image_arr[i][j]) + + return image_arr + + +"""color temperature table for white balancing""" +kelvin_table = { + 1000: (255,56,0), + 1500: (255,109,0), + 2000: (255,137,18), + 2500: (255,161,72), + 3000: (255,180,107), + 3500: (255,196,137), + 4000: (255,209,163), + 4500: (255,219,186), + 5000: (255,228,206), + 5500: (255,236,224), + 6000: (255,243,239), + 6500: (255,249,253), + 7000: (245,243,255), + 7500: (235,238,255), + 8000: (227,233,255), + 8500: (220,229,255), + 9000: (214,225,255), + 9500: (208,222,255), + 10000: (204,219,255)} + + +"""white balance""" +def convert_temp(image, temp): + r, g, b = kelvin_table[temp] + matrix = ( r / 255.0, 0.0, 0.0, 0.0, + 0.0, g / 255.0, 0.0, 0.0, + 0.0, 0.0, b / 255.0, 0.0 ) + return image.convert('RGB', matrix) + + +"""get the image output""" +def make_image(file_dir, filename): + + x, y, c = read_hit_file(file_dir) + + for i in range(len(c[0])): + + if c[0][i] > threshold_e: + c[0][i] = np.random.normal(c[0][i], smearing) + if c[0][i] < 0: + c[0][i] = 1 + else: + c[0][i] = 0 + + + image = np.zeros((1944, 2592)) + for i in range(len(x[0])): + image[len(image) - int(y[0][i]) - 1][int(x[0][i])] = c[0][i] + + image = image * conversion_factor + + maxY, maxX = get_brightest_pixel(image) + + if pass_edge_check(maxX, maxY, (2592, 1944)) == True: + + x0, x1, y0, y1 = get_crop_range(maxX, maxY) + + image = image[y0:y1, x0:x1] + + for i in range(len(image)): + for j in range(len(image[i])): + if image[i][j] < 0: + image[i][j] = 0 + + if image[i][j] > 255: + image[i][j] = 255.0 + + image = color_processing(image) + + image = Image.fromarray(np.uint8(np.array(image)), mode='RGB') + + image = convert_temp(image, 5000) + + image = np.array(image) + + image_gray = np.zeros((64, 64)) + + for i in range(len(image)): + for j in range(len(image[i])): + image_gray[i][j] = 0.299 * image[i][j][0] + 0.587 * image[i][j][1] + 0.114 * image[i][j][2] + + image = image_gray + + image = background_addition(image) + + if not os.path.exists("./output/" + file_qdc + "/" + str(particle_type) + "_2Dimages_fac" + str(conversion_factor) + "/"): + os.mkdir("./output/" + file_qdc + "/" + str(particle_type) + "_2Dimages_fac" + str(conversion_factor) + "/") + + f = open("./output/" + file_qdc + "/" + str(particle_type) + "_2Dimages_fac" + str(conversion_factor) + "/" + filename, 'w') + + for i in range(len(image)): + for k in range(len(image[i])): + f.write(str(int(image[i][k])) + " ") + f.write("\n") + f.close() + + +particle_type = "mu+" + +conversion_factor = 2 + +threshold_e = 0 + +smearing = 0 + +file_qdc = "qdc_res_0_smear_0_slope_0" +model = 'A510_ISO881' + + +for filename in os.listdir("./output/" + file_qdc + "/" + str(particle_type) + "/"): + + if os.path.exists("./output/" + file_qdc + "/" + str(particle_type) + "_2Dimages_fac" + str(conversion_factor) + "/" + filename): + pass + else: + print("working on " + str(filename)) + make_image("./output/" + file_qdc + "/" + str(particle_type) + "/" + filename, filename) + diff --git a/Simultaion_with_CNN/MC_simulation/Simulator_interface.py b/Simultaion_with_CNN/MC_simulation/Simulator_interface.py new file mode 100644 index 0000000..ab3efdf --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/Simulator_interface.py @@ -0,0 +1,165 @@ +r'''Class that handles various inputs and systematic parameters + and simulates muons''' +import os, sys +import subprocess +from pipes import quote +from math import * +import numpy as np + +class DECOSimulator(): + r'''Simulator class for muons that uses allpix + and GEANT''' + + def __init__(self, pid, energy, pos, theta, **kwargs): + self.pid = pid #Particle id, ie mu+, e- + self.pos = pos + self.energy = energy + self.theta = float(theta) + self.phi = float(kwargs.pop('phi', 0.)) + self.depletion_thickness = float(kwargs.pop('depletion_thickness', 26.3)) + #todo: NOTE: in unit of um !!!!!!! + # Set other possible systematics with kwargs, including + # electric fields, pixel size, etc. + + self.path_to_root = os.getenv('DECO_ROOT_PATH') + self.path_to_geant = os.getenv('DECO_GEANT_PATH') + self.path_to_allpix = os.getenv('DECO_ALLPIX_PATH') + self.base_command = self.path_to_allpix + ' -c ./htc_wildfire/source_measurement.conf -o' + + self.qdc_resolution = int(kwargs.pop('qdc_resolution', 0.)) + self.qdc_smearing = int(kwargs.pop('qdc_smearing', 0.)) + self.qdc_slope = kwargs.pop('qdc_slope', 0.) + + self.electronics_noise = kwargs.pop('electronics_noise', 0.) + self.bias_v = kwargs.pop('bias_v', -15) + self.gain = kwargs.pop('gain', 1.0) + + self.pixel_size = float(kwargs.pop('pixel_size', 1.0)) + self.seed = int(kwargs.pop('seed', -1)) + + def write_source_file(self, n_events, want_plot): + + with open('./htc_wildfire/source_measurement_replace.conf', 'r') as f: + data = f.readlines() + data[3] = data[3].format(n_events) + + if self.seed != -1: + data[3] = data[3] + "random_seed = " + str(self.seed) + "\n" + + data[18] = data[18].format(str(self.pos[0]) + " " + str(self.pos[1]) + " " + str(self.pos[2]) + "um") + + theta = radians(self.theta) + phi = radians(self.phi) + + print(self.theta, self.phi, self.energy) + dirx = -1*sin(theta)*cos(phi) + diry = -1*sin(theta)*sin(phi) + dirz = -1*cos(theta) + + data[20] = data[20].format(str(dirx) + " " + str(diry) + " " + str(dirz)) + + if want_plot == 'true': + data[38] = data[38].format('true' + "\n" + "output_linegraphs = true \n output_plots_step = 10ps \n output_plots_align_pixels = true \n output_plots_use_pixel_units = true") + else: + data[38] = data[38].format("false") + + data[50] = data[50].format(str(self.qdc_resolution)) + data[51] = data[51].format(str(self.qdc_smearing)) + data[52] = data[52].format(str(self.qdc_slope)) + data[53] = data[53].format(str(self.gain)) + data[54] = data[54].format(str(self.electronics_noise)) + data[28] = data[28].format(str(self.bias_v)) + data[29] = "\n" #data[29].format(str(self.depletion_thickness)) + + with open('./htc_wildfire/source_measurement.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + + def write_detector_file(self): + + with open('./htc_wildfire/htc_wildfire_shielded_replace.conf', 'r') as f: + data = f.readlines() + data[3] = 'pixel_size = ' + str(self.pixel_size) + "um " + str(self.pixel_size) + "um\n" + data[4] = data[4].format(self.depletion_thickness) + with open('./htc_wildfire/htc_wildfire_shielded.conf', 'w') as wf: + wf.writelines(data) + wf.close() + + + def set_output_file_name(self): + # write unique file name depending on parameters + if not os.path.exists("./output"): + os.system("mkdir ./output") + if not os.path.exists("./output/" + "qdc_res_" + str(self.qdc_resolution) + "_smear_" + + str(self.qdc_smearing) + "_slope_" + str(self.qdc_slope) + "/"): + os.system("mkdir ./output/" + "qdc_res_" + str(self.qdc_resolution) + "_smear_" + + str(self.qdc_smearing) + "_slope_" + str(self.qdc_slope) + "/") + if not os.path.exists("./output/" + "qdc_res_" + str(self.qdc_resolution) + "_smear_" + + str(self.qdc_smearing) + "_slope_" + str(self.qdc_slope) + "/" + str(self.pid)): + os.system("mkdir ./output/" + "qdc_res_" + str(self.qdc_resolution) + "_smear_" + + str(self.qdc_smearing) + "_slope_" + str(self.qdc_slope) + "/" + str(self.pid)) + + self.outfile = "qdc_res_" + str(self.qdc_resolution) + "_smear_" + str(self.qdc_smearing) + "_slope_" + str(self.qdc_slope) + "/" + str(self.pid) + \ + "/" + "{}_theta_{:.1f}_phi_{:.1f}_thickiness_{:.1f}_highstats.txt".format(self.energy, self.theta, self.phi, self.depletion_thickness) + pass + + + def get_output_file_name(self): + try: + return self.outfile + except: + self.set_output_file_name() + return self.outfile + + + def run_simulation(self, n_events, want_charge_plot='false'): + + self.source_local_env() + # Run the allpix simulation + + output_file = self.get_output_file_name() + self.write_source_file(n_events, want_charge_plot) + self.write_detector_file() + + my_command = self.base_command[:] + + my_command += ' DepositionGeant4.particle_type="{}"'.format(self.pid) + my_command += ' -o DepositionGeant4.source_energy="{}"'.format(self.energy) + my_command += ' -o TextWriter.file_name="' + output_file + '"' + + """check if simulated using file name""" + if self.check_if_simulated(output_file) is True: + print("has been simulated") + return + + with open(os.devnull, 'w') as devnull: + subprocess.call(my_command, shell=True, stdout=devnull) + + return + + + def check_if_simulated(self, filename): + # Check to see if simulation has already been run for + # this set of parameters + + if not os.path.exists("./output/" + filename): + return False + + return True + + + + def source_local_env(self): + geant = self.path_to_geant + "/bin/geant4.sh" + root = self.path_to_root + "/bin/thisroot.sh" + + pass + # Run source scripts for ROOT and GEANT if the + # simulation doesn't work + + + + + + diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/README.md b/Simultaion_with_CNN/MC_simulation/htc_wildfire/README.md new file mode 100644 index 0000000..3a6b976 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/README.md @@ -0,0 +1,8 @@ +# DECO Simulation with particle beams + +This directory is for simulating the detector described in the DECO measurement on the depletion thickness in cell phone camera image sensors. + +A Monolithic pixel is used, with a depletion thickness of 26.3um and extra non-instrumented silicon is added behind the chip to aid in the production of worm like signatures in the detector. +No misalignment is added but the absolute position and orientation of the detector is specified. + +The setup of the simulation chain includes Generic Digitization, instantiation of a particle beam, as well as detector geometry relative to the incident beam. There is a python script which can be used to submit batches of beams to cover large phase spaces with ease. diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/allpix_example_detector.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/allpix_example_detector.conf new file mode 100644 index 0000000..07eb855 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/allpix_example_detector.conf @@ -0,0 +1,31 @@ +# Reference http://www.sciencedirect.com/science/article/pii/S0168900210012982 +type = "hybrid" + +number_of_pixels = 256 256 +pixel_size = 55um 55um + +sensor_thickness = 300um +sensor_excess = 1mm + +bump_sphere_radius = 9.0um +bump_cylinder_radius = 7.0um +bump_height = 20.0um + +chip_thickness = 100um +chip_excess_top = 1610um +chip_excess_bottom = 1610um +chip_excess_right = 10um +chip_excess_left = 10um + +[support] +thickness = 1.76mm +size = 47mm 79mm +offset = 0 -22.25mm +material = "g10" + +[support] +thickness = 0.8mm +size = 15mm 15mm +material = "aluminum" +location = "absolute" +offset = 0 0 -3mm diff --git a/htc_wildfire/detector_replace.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/detector.conf similarity index 80% rename from htc_wildfire/detector_replace.conf rename to Simultaion_with_CNN/MC_simulation/htc_wildfire/detector.conf index e95c9f8..8f9cdec 100644 --- a/htc_wildfire/detector_replace.conf +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/detector.conf @@ -2,5 +2,5 @@ type = "htc_wildfire_shielded" position = 0um 0um 0um orientation_mode = "zyx" -orientation = {}deg {}deg {}deg +orientation = 0.0deg 0deg 0deg #orientation is azimuth then zenith angle diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded.conf new file mode 100644 index 0000000..e00c7db --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded.conf @@ -0,0 +1,16 @@ +type = "monolithic" + +number_of_pixels = 2592 1944 +pixel_size = 0.9um 0.9um +sensor_thickness = 26.3um +#sensor_excess_direction = 'bottom' +#sensor_excess = 0.5mm +chip_thickness = 1mm + +#For shielding +# [support] +# thickness = 1.0mm +# size = 35mm 35mm +# material = "lead" +# location = "absolute" +# offset = 0 0 -20mm diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded_replace.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded_replace.conf new file mode 100644 index 0000000..bcd42d4 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/htc_wildfire_shielded_replace.conf @@ -0,0 +1,16 @@ +type = "monolithic" + +number_of_pixels = 2592 1944 +pixel_size = {}um {}um +sensor_thickness = {}um +#sensor_excess_direction = 'bottom' +#sensor_excess = 0.5mm +chip_thickness = 1mm + +#For shielding +# [support] +# thickness = 1.0mm +# size = 35mm 35mm +# material = "lead" +# location = "absolute" +# offset = 0 0 -20mm diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/readme.txt b/Simultaion_with_CNN/MC_simulation/htc_wildfire/readme.txt new file mode 100644 index 0000000..042801e --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/readme.txt @@ -0,0 +1 @@ +These are the Allpix2 and GEANT4 setups for this simuation diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement.conf new file mode 100644 index 0000000..3cb90d9 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement.conf @@ -0,0 +1,77 @@ +[AllPix] +log_level = "WARNING" +log_format = "DEFAULT" +number_of_events = 1 +detectors_file = "detector.conf" +model_paths = ./ + +# Set surrounding material to air +[GeometryBuilderGeant4] +world_material = "vacuum" + +# Beam parameters and typical GEANT4 physics lists +[DepositionGeant4] +physics_list = FTFP_BERT_LIV +#enable_pai = true +charge_creation_energy = 3.62eV +#particle_type = "mu+" +#source_energy = 1GeV +source_position = 0 0 13.15um +source_type = "beam" +beam_direction = -0.612625791461 0.404695865787 -0.678904187535 +beam_divergence = 0.1mrad 0.1mrad +max_step_length = 1um +output_plots = false + +# Optional electric field applied across pixels +[ElectricFieldReader] +model="linear" +bias_voltage=-20V +depletion_voltage=-10V + +# Pick between this and GenericPropagation (do holes or e- move)? +[GenericPropagation] +temperature = 293K +propagate_electrons = true +propagate_holes = true +charge_per_step = 1 +integration_time = 50ns +output_plots = false + +# combines single charges to sets of charges for processing +[SimpleTransfer] +max_depth_distance = 1.0um + +# Add noise to the other pixels, set thresholds for digitization +# Smear signal, add gain +[DefaultDigitizer] +output_plots = false +threshold = 1e +threshold_smearing = 1e +qdc_resolution = 0 +qdc_smearing = 0e +qdc_slope = 7e +gain = 1.0 +electronics_noise = 5e +allow_zero_qdc = false + +# What to write to output +# [ROOTObjectWriter] +# exclude = DepositedCharge, PropagatedCharge +# file_name = "output_measurement_3eV.root" + +[TextWriter] +include = "PixelHit" + +# Visualization stuff +#[VisualizationGeant4] +#mode = "gui" +#view_style = "surface" +###simple_view = false +#transparency = 0.3 +#trajectories_color_mode = "charge" + +# The object writer listens to all output data +#[ROOTObjectWriter] +# specify the output file (default file name is used if omitted) +#file_name = "data.root" diff --git a/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement_replace.conf b/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement_replace.conf new file mode 100644 index 0000000..e49a3e7 --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/htc_wildfire/source_measurement_replace.conf @@ -0,0 +1,77 @@ +[AllPix] +log_level = "WARNING" +log_format = "DEFAULT" +number_of_events = {} +detectors_file = "detector.conf" +model_paths = ./ + +# Set surrounding material to air +[GeometryBuilderGeant4] +world_material = "vacuum" + +# Beam parameters and typical GEANT4 physics lists +[DepositionGeant4] +physics_list = FTFP_BERT_LIV +#enable_pai = true +charge_creation_energy = 3.62eV +#particle_type = "mu+" +#source_energy = 1GeV +source_position = {} +source_type = "beam" +beam_direction = {} +beam_divergence = 0.1mrad 0.1mrad +max_step_length = 1um +output_plots = false + +# Optional electric field applied across pixels +[ElectricFieldReader] +model="linear" +bias_voltage={}V +depletion_thickness={}um + +# Pick between this and GenericPropagation (do holes or e- move)? +[GenericPropagation] +temperature = 293K +propagate_electrons = true +propagate_holes = false +charge_per_step = 1 +integration_time = 50ns +output_plots = {} + +# combines single charges to sets of charges for processing +[SimpleTransfer] +max_depth_distance = 0.5um + +# Add noise to the other pixels, set thresholds for digitization +# Smear signal, add gain +[DefaultDigitizer] +output_plots = false +threshold = 1e +threshold_smearing = 0e +qdc_resolution = {} +qdc_smearing = {}e +qdc_slope = {}e +gain = {} +electronics_noise = {}e +allow_zero_qdc = false + +# What to write to output +# [ROOTObjectWriter] +# exclude = DepositedCharge, PropagatedCharge +# file_name = "output_measurement_3eV.root" + +[TextWriter] +include = "PixelHit" #, "MCTrack", "MCParticle" + +# Visualization stuff +#[VisualizationGeant4] +#mode = "gui" +#view_style = "surface" +###simple_view = false +#transparency = 0.3 +#trajectories_color_mode = "charge" + +# The object writer listens to all output data +#[ROOTObjectWriter] +# specify the output file (default file name is used if omitted) +#file_name = "data.root" diff --git a/Simultaion_with_CNN/MC_simulation/random_e_gamma_generator.py b/Simultaion_with_CNN/MC_simulation/random_e_gamma_generator.py new file mode 100644 index 0000000..716c76b --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/random_e_gamma_generator.py @@ -0,0 +1,35 @@ +#!/usr/bin/env python + +import numpy as np +from math import * + + +"""Weight function corresponds to isotropic zenith angle distribution with detector effective area effect""" +def checker(ang): + + checker = np.random.random() + + val = cos(ang) * (1 - cos(ang))/2 + + if checker <= val: + return True + + else: + return False + + +def get_rand_ini_e_gamma(): + + while 1: # and not zenith_dir > pi/4: + + zen_ang = np.random.random() * pi / 2 + + if checker(zen_ang) == True: + + """---------azimuth part---------""" + azimuth_dir = np.random.random() * 2 * pi + + return zen_ang, azimuth_dir + + + diff --git a/Simultaion_with_CNN/MC_simulation/random_muon_generator.py b/Simultaion_with_CNN/MC_simulation/random_muon_generator.py new file mode 100644 index 0000000..d76684e --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/random_muon_generator.py @@ -0,0 +1,44 @@ +#!/usr/bin/env python + +import numpy as np +from math import * + + +"""Weight function is according to the atmospheric muon flux and detector effective area""" +def check_in_range(E, theta, rand): + val = pow(E, -2.7) * ( 1/ (1 + 1.11 * E * np.cos(theta) /115) + 0.054 / ( 1 + 1.11 * E * np.cos(theta) /850) ) * pow(E / (E + 2/cos(theta)), 2.7) * cos(theta) * sin(theta) + + if rand < val: + return True + else: + return False + + +def get_rand_ini_muon(E_min, E_max): + E_min = E_min + 0.105658 + E_max = E_max + 0.105658 + + factor = 1.1 * pow(E_min, -2.7) + + checker = False + + + zenith_dir = 0 + energy = 0 + + while not checker: # and not zenith_dir > pi/4: + energy = np.random.random() * (E_max - E_min) + E_min + zenith_dir = np.random.random() * (np.radians(90)) + curr = np.random.random() * factor + + checker = check_in_range(energy, zenith_dir, curr) + + """---------azimuth part---------""" + azimuth_dir = np.random.random() * 2 * pi + + + + return energy - 0.105658, zenith_dir, azimuth_dir + + + diff --git a/Simultaion_with_CNN/MC_simulation/readme.txt b/Simultaion_with_CNN/MC_simulation/readme.txt new file mode 100644 index 0000000..cb22a4b --- /dev/null +++ b/Simultaion_with_CNN/MC_simulation/readme.txt @@ -0,0 +1,13 @@ +Run with: + +source geant/bin/geant.sh +source root/bin/thisroot.sh + +python2 deco_rand_muon_simulation.py + +The files generated are charge deposited on pixels, run Single_plot.py to convert them into luminance on image +Note that background files extracted from real events should be put into appropriate folder or set to 0 + +The folder htcwildfire contains all parameters for the HTC wildfire model. In the file source_measurement_replace.conf +"MCTrack", "MCParticle" has been commented out, and they represent the MC true value of particles produced +in this simulation (remove the '#' to make them work) diff --git a/Simultaion_with_CNN/Plot_parameters/plot_hillas.py b/Simultaion_with_CNN/Plot_parameters/plot_hillas.py new file mode 100644 index 0000000..170cd88 --- /dev/null +++ b/Simultaion_with_CNN/Plot_parameters/plot_hillas.py @@ -0,0 +1,195 @@ +import numpy as np +import os, sys +import matplotlib.pylab as plt +plt.style.use('IceCube') +import matplotlib.pyplot as plt +from math import * + +# from Brent, calculate Hillas parameters of events +#charge_coords = [[x_coords], [y_coords], [charges]] +def hillas(charge_coords): + #print(charge_coords.shape) + x = 0 + y = 0 + x2 = 0 + y2 = 0 + xy = 0 + CHARGE = 0 + #print(charge_coords.shape) + CHARGE = np.sum(charge_coords[2]) + x = np.sum(charge_coords[0] * charge_coords[2]) + y = np.sum(charge_coords[1] * charge_coords[2]) + x2 = np.sum(charge_coords[0] ** 2 * charge_coords[2]) + y2 = np.sum(charge_coords[1] ** 2 * charge_coords[2]) + xy = np.sum(charge_coords[0] * charge_coords[1] * charge_coords[2]) + x /= CHARGE + y /= CHARGE + x2 /= CHARGE + y2 /= CHARGE + xy /= CHARGE + S2_x = x2 - x ** 2 + S2_y = y2 - y ** 2 + S_xy = xy - x * y + d = S2_y - S2_x + a = (d + np.sqrt(d ** 2 + 4 * S_xy ** 2)) / (2 * S_xy) + b = y - a * x + width = np.sqrt((S2_y + a ** 2 * S2_x - 2 * a * S_xy) / (1 + a ** 2)) + length = np.sqrt((S2_x + a ** 2 * S2_y + 2 * a * S_xy) / (1 + a ** 2)) + miss = np.abs(b / np.sqrt(1 + a ** 2)) + dis = np.sqrt(x ** 2 + y ** 2) + q_coord = (x - charge_coords[0]) * (x / dis) + (y - charge_coords[1]) * (y / dis) + q = np.sum(q_coord * charge_coords[2]) / CHARGE + q2 = np.sum(q_coord ** 2 * charge_coords[2]) / CHARGE + azwidth = q2 - q ** 2 + + return [width, length, miss, dis, azwidth] + + +# convert a 64x64 charge deposited plot or real deco image into shape of [[x_coords], [y_coords], [charges]] +def convert_hillas(filename): + f = open(filename, 'r') + + y_list = [] + x_list = [] + charge_list = [] + + y_pos = 0.0 + + while 1: + + temp = f.readline().split() + + if len(temp) < 1: + break + + x_pos = 0.0 + + for num in temp: + x_list.append(x_pos) + y_list.append(y_pos) + + if float(num) >= 5.0: + charge_list.append(float(num)) + else: + charge_list.append(0.0) + + x_pos += 1.0 + + y_pos += 1.0 + + y_list = len(y_list) - np.array(y_list) - 1 + x_list = np.array(x_list) + charge_list = np.array(charge_list) + + + return [x_list, y_list, charge_list] + + +#print hillas(convert_hillas("./output/mu+_2Dimages/1.0004537985GeV_theta_57.0_phi_75.9_thickiness_26.3_highstats.txt")) + +""" +extracting cnn classified data +""" +width_list_cnn = [] +length_list_cnn = [] +count = 0 + + +for filename in os.listdir("../MC_z_clean_img/real_clean_image_thres10/"): + if filename.endswith(".txt"): + print(count) + count += 1 + try : + result = hillas(convert_hillas("../MC_z_clean_img/real_clean_image_thres10/" + filename)) + + width_list_cnn.append(log10(result[0])) + length_list_cnn.append(log10(result[1])) + # if log10(result[0]) > 0.4: + # print filename + except: + pass + + + + +""" +extracting simulated data +""" + +def getdata_sim(file_dir): + width_list_sim3 = [] + length_list_sim3 = [] + count = 0 + + for filename in os.listdir(file_dir): + + if filename.endswith(".txt"): + + try: + result = hillas(convert_hillas(file_dir + "/" + filename)) + + width_list_sim3.append(log10(result[0])) + length_list_sim3.append(log10(result[1])) + + count += 1 + except: + pass + + return width_list_sim3, length_list_sim3, count + + + +folder_dir = "qdc_res_0_smear_0_slope_4" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +width_list_sim5, length_list_sim5, count5 = getdata_sim(file_dir) + + +folder_dir = "qdc_res_0_smear_0_slope_3" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +width_list_sim25, length_list_sim25, count25 = getdata_sim(file_dir) + + +folder_dir = "qdc_res_0_smear_0_slope_5" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +width_list_sim100, length_list_sim100, count100 = getdata_sim(file_dir) + + +f = plt.figure() + +plt.hist(width_list_cnn, bins=np.linspace(-0.2, 0.7, 41), weights=np.zeros(len(width_list_cnn)) + (1 / count), label='real events', alpha=0.5) + +plt.hist(width_list_sim5, bins=np.linspace(-0.2, 0.7, 41), weights=np.zeros(len(width_list_sim5)) + (1 / count5), label='simulation with bias voltage -5V', histtype='step', lw=5) +plt.hist(width_list_sim25, bins=np.linspace(-0.2, 0.7, 41), weights=np.zeros(len(width_list_sim25)) + (1 / count25), label='simulation with bias voltage -25V', histtype='step', lw=5, color='black') +plt.hist(width_list_sim100, bins=np.linspace(-0.2, 0.7, 41), weights=np.zeros(len(width_list_sim100)) + (1 / count100), label='simulation with bias voltage -100V', histtype='step', lw=5) + + +plt.legend() +plt.xlabel("log10(width / pixel)") +plt.ylabel("# images normalized") +#plt.xscale('log') +#plt.yscale('log') + + + +f2 = plt.figure() + +plt.hist(length_list_cnn, bins=np.linspace(0, 1.5, 31), weights=np.zeros(len(width_list_cnn)) + (1 / count), label='real events', alpha=0.5) + +plt.hist(length_list_sim5, bins=np.linspace(0, 1.5, 31), weights=np.zeros(len(width_list_sim5)) + (1 / count5), label='simulation with bias voltage -5V', histtype='step', lw=5) +plt.hist(length_list_sim25, bins=np.linspace(0, 1.5, 31), weights=np.zeros(len(width_list_sim25)) + (1 / count25), label='simulation with bias voltage -25V', histtype='step', lw=5, color='black') +plt.hist(length_list_sim100, bins=np.linspace(0, 1.5, 31), weights=np.zeros(len(width_list_sim100)) + (1 / count100), label='simulation with bias voltage -100V', histtype='step', lw=5) + + +plt.legend() +plt.xlabel("log10(length / pixel)") +plt.ylabel("# images normalized") +#plt.yscale('log') +#plt.xscale('log') + + + +plt.show() + diff --git a/Simultaion_with_CNN/Plot_parameters/plot_moments.py b/Simultaion_with_CNN/Plot_parameters/plot_moments.py new file mode 100644 index 0000000..168d71b --- /dev/null +++ b/Simultaion_with_CNN/Plot_parameters/plot_moments.py @@ -0,0 +1,101 @@ +import matplotlib.pylab as plt +from math import * +import numpy as np +plt.style.use("IceCube") +import os, sys + +def getdata(filename): + f = open(filename, 'r') + + output = [] + + while 1: + temp = f.readline().split() + + if len(temp) < 1: + break + + for num in temp: + if float(num) != 0: + output.append(float(num)) + + return np.array(output) + + + + +count = 0 +x_sig_real = [] +x_mean_real = [] +for filename in os.listdir("../MC_z_clean_img/real_clean_image_thres10/"): + output = getdata("../MC_z_clean_img/real_clean_image_thres10/" + filename) + x_sig_real.append(np.sqrt(np.power(output, 2).mean() - output.mean()**2)) + x_mean_real.append(output.mean()) + count += 1 + print(count) + + +def getdata_sim(file_dir): + + count = 0 + x_sig_sim = [] + x_mean_sim = [] + + for filename in os.listdir(file_dir): + output = getdata(file_dir + "/" + filename) + x_sig_sim.append(np.sqrt(np.power(output, 2).mean() - output.mean()**2)) + x_mean_sim.append(output.mean()) + count += 1 + print(count) + + return x_sig_sim, x_mean_sim, count + + +folder_dir = "qdc_res_0_smear_0_slope_4" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +x_sig_sim5, x_mean_sim5, count5 = getdata_sim(file_dir) + + +folder_dir = "qdc_res_0_smear_0_slope_3" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +x_sig_sim25, x_mean_sim25, count25 = getdata_sim(file_dir) + + +folder_dir = "qdc_res_0_smear_0_slope_5" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +x_sig_sim100, x_mean_sim100, count100 = getdata_sim(file_dir) + + +f = plt.figure() + +plt.hist(x_sig_real, bins=np.linspace(0, 80, 41), weights=np.zeros(len(x_sig_real)) + (1 / count), label='real events', alpha=0.5) + +plt.hist(x_sig_sim5, bins=np.linspace(0, 80, 41), weights=np.zeros(len(x_sig_sim5)) + (1 / count5), label='simulation with bias voltage -5V', histtype='step', lw=5) +plt.hist(x_sig_sim25, bins=np.linspace(0, 80, 41), weights=np.zeros(len(x_sig_sim25)) + (1 / count25), label='simulation with bias voltage -25V', histtype='step', lw=5, color='black') +plt.hist(x_sig_sim100, bins=np.linspace(0, 80, 41), weights=np.zeros(len(x_sig_sim100)) + (1 / count100), label='simulation with bias voltage -100V', histtype='step', lw=5) + +plt.legend() +plt.xlabel(r"$\sigma_L$") +plt.ylabel("# images normalized") + + + + +f2 = plt.figure() + +plt.hist(np.array(x_sig_real)/np.array(x_mean_real), bins=np.linspace(0.5, 1.5, 41), weights=np.zeros(len(x_sig_real)) + (1 / count), label='real events', alpha=0.5) + +plt.hist(np.array(x_sig_sim5)/np.array(x_mean_sim5), bins=np.linspace(0.5, 1.5, 41), weights=np.zeros(len(x_sig_sim5)) + (1 / count5), label='simulation with bias voltage -5V', histtype='step', lw=5) +plt.hist(np.array(x_sig_sim25)/np.array(x_mean_sim25), bins=np.linspace(0.5, 1.5, 41), weights=np.zeros(len(x_sig_sim25)) + (1 / count25), label='simulation with bias voltage -25V', histtype='step', lw=5, color='black') +plt.hist(np.array(x_sig_sim100)/np.array(x_mean_sim100), bins=np.linspace(0.5, 1.5, 41), weights=np.zeros(len(x_sig_sim100)) + (1 / count100), label='simulation with bias voltage -100V', histtype='step', lw=5) + +plt.legend() +plt.xlabel(r"$\frac{\sigma_L}{}$") +plt.ylabel("# images normalized") + +plt.show() + + diff --git a/Simultaion_with_CNN/Plot_parameters/plot_pixel_hist.py b/Simultaion_with_CNN/Plot_parameters/plot_pixel_hist.py new file mode 100644 index 0000000..87c3ca1 --- /dev/null +++ b/Simultaion_with_CNN/Plot_parameters/plot_pixel_hist.py @@ -0,0 +1,141 @@ +import matplotlib.pylab as plt +from math import * +import numpy as np +plt.style.use("IceCube") +import os, sys + +def getdata(filename, if_sim): + f = open(filename, 'r') + + output = [] + + while 1: + temp = f.readline().split() + + if len(temp) < 1: + break + + for num in temp: + if float(num) > 0: + output.append(float(num)) + + return np.array(output), np.array(output).max() + + +output_total_cnn = np.array([]) +output_max_cnn = [] +output_sum = [] +count_real = 0 + +for filename in os.listdir("../MC_z_clean_img/real_clean_image_thres10/"): + if filename.endswith(".txt"): + + output, max_lum = getdata("../MC_z_clean_img/real_clean_image_thres10/" + filename, False) + output_max_cnn.append(max_lum) + output_sum.append(output.sum()) + if output.sum() > 10000: + print(filename) + output_total_cnn = np.concatenate((output_total_cnn, output)) + count_real += 1 + +def getdata_sim(file_dir): + + output_total_sim2 = np.array([]) + output_max_sim2 = [] + output_sum_sim2 = [] + + count_1 = 0 + + for filename in os.listdir(file_dir): + if filename.endswith(".txt"): + + print(count_1) + + count_1 += 1 + + + output, max_lum = getdata(file_dir + "/" + filename, True) + output_max_sim2.append(max_lum) + output_sum_sim2.append(output.sum()) + + output_total_sim2 = np.concatenate((output_total_sim2, output)) + + return output_total_sim2, output_max_sim2 ,output_sum_sim2, count_1 + + +folder_dir = "qdc_res_0_smear_0_slope_4" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +output_total_sim5, output_max_sim5 ,output_sum_sim5, count5 = getdata_sim(file_dir) + + + +folder_dir = "qdc_res_0_smear_0_slope_3" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +output_total_sim25, output_max_sim25 ,output_sum_sim25, count25 = getdata_sim(file_dir) + + +folder_dir = "qdc_res_0_smear_0_slope_5" +source_dir = "mu+_2Dimages_fac2" +file_dir = "../MC_z_clean_img/sim_clean_image_thres10/" + str(folder_dir) + "/" + str(source_dir) +output_total_sim100, output_max_sim100 ,output_sum_sim100, count100 = getdata_sim(file_dir) + + + + +f = plt.figure() + +plt.hist(output_total_cnn, bins=np.linspace(0, 255, 86), weights=np.zeros(len(output_total_cnn)) + (1 / count_real), alpha=0.5, label='real events') + +plt.hist(output_total_sim5, bins=np.linspace(0, 255, 86), weights=np.zeros(len(output_total_sim5)) + (1 / count5), histtype='step', label='simulation with bias voltage -5V', lw=5) +plt.hist(output_total_sim25, bins=np.linspace(0, 255, 86), weights=np.zeros(len(output_total_sim25)) + (1 / count25), histtype='step', label='simulation with bias voltage -25V', lw=5, color='black') +plt.hist(output_total_sim100, bins=np.linspace(0, 255, 86), weights=np.zeros(len(output_total_sim100)) + (1 / count100), histtype='step', label='simulation with bias voltage -100V', lw=5) + +plt.legend() +plt.xlabel("Luminance on pixel") +plt.yscale('log') +plt.ylabel("# pixels normalized") + + + + +f2 = plt.figure() + +plt.hist(output_max_cnn, bins=np.linspace(0, 255, 26), weights=np.zeros(len(output_max_cnn)) + (1 / count_real), alpha=0.5, label='real events') + +plt.hist(output_max_sim5, bins=np.linspace(0, 255, 26), weights=np.zeros(len(output_max_sim5)) + (1 / count5), histtype='step', label='simulation with bias voltage -5V', lw=5) +plt.hist(output_max_sim25, bins=np.linspace(0, 255, 26), weights=np.zeros(len(output_max_sim25)) + (1 / count25), histtype='step', label='simulation with bias voltage -25V', lw=5, color='black') +plt.hist(output_max_sim100, bins=np.linspace(0, 255, 26), weights=np.zeros(len(output_max_sim100)) + (1 / count100), histtype='step', label='simulation with bias voltage -100V', lw=5) + +plt.legend() +plt.xlabel("Brightest luminance on pixel") +plt.ylabel("# images normalized") +#plt.yscale('log') + + + + + + +f3 = plt.figure() + + + +plt.hist(output_sum, bins=np.logspace(3, 4.5, 41), weights=np.zeros(len(output_sum)) + (1 / count_real), alpha=0.5, label='real events') + +plt.hist(output_sum_sim5, bins=np.logspace(3, 4.5, 41), weights=np.zeros(len(output_sum_sim5)) + (1 / count5), histtype='step', label='simulation with bias voltage -5V', lw=5) +plt.hist(output_sum_sim25, bins=np.logspace(3, 4.5, 41), weights=np.zeros(len(output_sum_sim25)) + (1 / count25), histtype='step', label='simulation with bias voltage -25V', lw=5, color='black') +plt.hist(output_sum_sim100, bins=np.logspace(3, 4.5, 41), weights=np.zeros(len(output_sum_sim100)) + (1 / count100), histtype='step', label='simulation with bias voltage -100V', lw=5) + +plt.xscale('log') +plt.legend() +plt.xlabel("Luminance summation") +plt.ylabel("# images normalized") + + + + + +plt.show() + diff --git a/Simultaion_with_CNN/Plot_parameters/readme.txt b/Simultaion_with_CNN/Plot_parameters/readme.txt new file mode 100644 index 0000000..7d9ab0e --- /dev/null +++ b/Simultaion_with_CNN/Plot_parameters/readme.txt @@ -0,0 +1 @@ +This folder contains files to plot the Hillas parameters, Luminance distirbution, and moments of luminance distribution diff --git a/Simultaion_with_CNN/image_cleaning.py b/Simultaion_with_CNN/image_cleaning.py new file mode 100644 index 0000000..5ae5d5c --- /dev/null +++ b/Simultaion_with_CNN/image_cleaning.py @@ -0,0 +1,89 @@ +import numpy as np +from math import * +import matplotlib.pylab as plt +from PIL import Image +import os, sys + + +def getImage(filename): + + f_reader = open(filename, 'r') + image = [] + while 1: + temp = f_reader.readline().split() + temp_hold = [] + if len(temp) < 1: + break + + for content in temp: + temp_hold.append(float(content)) + + image.append(temp_hold) + + image = np.array(image) + + return image + + +def clean_img(filename): + + image = getImage(filename) + + noise_matrix = np.zeros((5, 5)) + threshold + + a = 0.2138 + b = 0.479 + c = 0.985 + aperture = [[0., a, b, a, 0.], + [a, c, 1., c, a], + [b, 1., 1., 1., b], + [a, c, 1., c, a], + [0., a, b, a, 0.]] + + aperture = np.array(aperture) + + image_pad = np.zeros((len(image) + 4, len(image[0]) + 4)) + + for i in range(len(image_pad))[2:-2]: + for j in range(len(image_pad[0]))[2:-2]: + + image_pad[i][j] = image[i-2][j-2] + + + for i in range(len(image_pad))[2:-2]: + for j in range(len(image_pad[0]))[2:-2]: + + luminance = np.sum(image_pad[i-2:i + 3, j-2:j + 3] * aperture) + + if luminance <= np.sum(noise_matrix * aperture): + image[i-2][j-2] = 0 + + + return image + + +# threshold of image cleaning +threshold = 10 + + +if not os.path.exists("./real_clean_worm_thres" + str(threshold)): + os.mkdir("./real_clean_worm_thres" + str(threshold)) + +file_dir = "../CNN_on_training_data/Worm_data/HTC_A510c/" + +for filename in os.listdir(file_dir): + + print(filename) + + clean_image = clean_img(file_dir + filename) + + f = open("./real_clean_worm_thres" + str(threshold) + "/" + filename, 'w') + + for i in range(len(clean_image)): + for j in range(len(clean_image[i])): + f.write(str(int(clean_image[i][j])) + " ") + + f.write("\n") + + + diff --git a/Simultaion_with_CNN/readme.txt b/Simultaion_with_CNN/readme.txt new file mode 100644 index 0000000..31158f8 --- /dev/null +++ b/Simultaion_with_CNN/readme.txt @@ -0,0 +1,5 @@ +This folder contains the Monte Carlo simulation of cosmic-ray muons and electrons/photons from radioactive decay in background detected by cell phone CMOS in MC_Simulation + +The image cleaning algorithm is included + +The CNN used to classify simulated events is in CNN_on_simulation diff --git a/Single_plot.py b/Single_plot.py new file mode 100644 index 0000000..0bc87ab --- /dev/null +++ b/Single_plot.py @@ -0,0 +1,156 @@ +import numpy as np +import pandas as pd +import matplotlib.pylab as plt +from matplotlib.colors import ListedColormap +import seaborn as sns +import os +plt.style.use("IceCube") +import matplotlib as mpl + +class DECOLeptonAnalyzer(): + r'''This class is for making plots like the ones in the + notebook that read in a bunch of simulated files + and make analysis level plots''' + + def __init__(self, pid, phi, thickness): + self.pid = pid + + self.thichness = thickness + + self.phi = phi + + + def read_hit_file(self, filename): + f = open(filename, 'r') + + xhits, yhits, charge = [], [], [] + + # skip first 2 lines + f.readline() + f.readline() + + x, y, c = [], [], [] + + flag = 1 + while 1: + + temp = f.readline().split() + + if len(temp) < 1 or temp[0] == '#': + break + + if temp[0] == '===': + continue + if temp[0] == '---': + if flag == 1: + flag = 0 + else: + xhits.append(x) + yhits.append(y) + charge.append(c) + x, y, c = [], [], [] + + else: + x.append(float(temp[1][:-1])) + y.append(float(temp[2][:-1])) + c.append(float(temp[3][:-1])) + + if len(x) > 0: + xhits.append(x) + yhits.append(y) + charge.append(c) + + return xhits, yhits, charge + + + + + def plot_single(self, en, numerica_e, ang, num_to_do): + x, y, c = self.read_hit_file( + "./output/{}/{}_theta_{}_phi_{}_thickiness_{}_highstats.txt".format(self.pid, en, float(ang), + float(self.phi), self.thichness)) + + for j in range(len(x))[:num_to_do]: + + image = np.zeros((4500, 4500)) + for i in range(len(x[j])): + image[int(y[j][i]), int(x[j][i])] = c[j][i] + + + med_x = np.median(x[j]) + med_y = np.median(y[j]) + size = 50. + + fig1 = plt.figure(1, figsize=(8, 8)) + ax = fig1.add_subplot(111) + fig1.set_facecolor('white') + + # my_cmap = ListedColormap(sns.color_palette("Blues", 50)) + # my_cmap = ListedColormap(sns.palplot(sns.cubehelix_palette(8, start=2, rot=0, dark=0, light=.95, reverse=True))) + my_cmap = mpl.cm.hot + + # image = np.where(image == 0.0, np.nan, image) + + im = ax.imshow(image, cmap=my_cmap, # interpolation="gaussian", + aspect="auto", vmax=100., vmin=0.0) + ax.set_xlim([med_x - size, med_x + size]) + ax.set_ylim([med_y - size, med_y + size]) + ax.set_xlabel("X (pixels)") + ax.set_ylabel("Y (pixels)") + + ax.grid(color="#ffffff") + cb = fig1.colorbar(im, orientation="vertical", + shrink=0.8, + fraction=0.05, + pad=0.15) + label = "Pixel Luminance" + cb.set_label(label) + ax.text(med_x + size * 0.3, med_y + size * 0.7, + "Simulation", fontsize=24, color='w', weight='heavy') + + if not os.path.exists("./2Dimages"): + os.mkdir("./2Dimages") + if not os.path.exists("./2Dimages/" + str(particle_type)): + os.mkdir("./2Dimages/" + str(particle_type)) + if not os.path.exists("./2Dimages/" + str(particle_type) + "/" + str(round(numerica_e/pow(10, 6), 4)) + "MeV"): + os.mkdir("./2Dimages/" + str(particle_type) + "/" + str(round(numerica_e/pow(10, 6), 4)) + "MeV") + + title = str(self.pid) + " initial energy is " + str(round(numerica_e/pow(10, 6), 4)) + "MeV\ntotal charge deposited is "\ + + str(round(np.array(c[j]).sum(), 3)) + "\ntotal deposited energy is: "\ + + str(round(np.array(c[j]).sum() * 3.62/pow(10, 6), 3)) + "MeV" + + ax.set_title(title) + + plt.savefig("./2Dimages/" + str(self.pid) + "/" + str(round(numerica_e/pow(10, 6), 4)) + "MeV/" + str(j) + ".png", bbox_inches='tight') + + plt.close() + + + + +energy = np.logspace(4, 10, 100) + +energy_list = [] + +for ene in energy: + energy_list.append(str(ene/pow(10, 6))+"MeV") + +theta = '45' +phi = '0' +particle_type = 'e-' + +thickness = 26.3 + +a = DECOLeptonAnalyzer(particle_type, phi, thickness) + +E_calibrate_list = [] +err_list = [] +E_exp_list = [] + +for i in range(len(energy)): + + print("working on " + str(round(energy[i]/pow(10, 6), 4)) + "MeV") + + a.plot_single(energy_list[i], energy[i], theta, 20) + + diff --git a/htc_wildfire/.DS_Store b/htc_wildfire/.DS_Store new file mode 100644 index 0000000..56e7c5d Binary files /dev/null and b/htc_wildfire/.DS_Store differ diff --git a/htc_wildfire/detector.conf b/htc_wildfire/detector.conf index 67b35da..8f9cdec 100644 --- a/htc_wildfire/detector.conf +++ b/htc_wildfire/detector.conf @@ -2,5 +2,5 @@ type = "htc_wildfire_shielded" position = 0um 0um 0um orientation_mode = "zyx" -orientation = 0.0deg 30.0deg 0deg +orientation = 0.0deg 0deg 0deg #orientation is azimuth then zenith angle diff --git a/htc_wildfire/htc_wildfire_shielded.conf b/htc_wildfire/htc_wildfire_shielded.conf index 606e11c..e00c7db 100644 --- a/htc_wildfire/htc_wildfire_shielded.conf +++ b/htc_wildfire/htc_wildfire_shielded.conf @@ -3,6 +3,7 @@ type = "monolithic" number_of_pixels = 2592 1944 pixel_size = 0.9um 0.9um sensor_thickness = 26.3um +#sensor_excess_direction = 'bottom' #sensor_excess = 0.5mm chip_thickness = 1mm @@ -13,84 +14,3 @@ chip_thickness = 1mm # material = "lead" # location = "absolute" # offset = 0 0 -20mm - -# Build the phone layout -[support] -# front screen -thickness = 3.0mm -size = 70mm 150mm -material = "plexiglass" -location = "absolute" -offset = 0 -50mm +3.5mm - -[support] -# lens -thickness = 4.0mm -size = 5mm 5mm -material = "plexiglass" -location = "absolute" -offset = 0 0 -3.5mm - -[support] -# back below lens -thickness = 3.0mm -size = 70mm 122.5mm -material = "aluminum" -location = "absolute" -offset = 0 -63.75mm -3.5mm - -[support] -# back above lens -thickness = 3.0mm -size = 70mm 22.5mm -material = "aluminum" -location = "absolute" -offset = 0 13.75mm -3.5mm - -[support] -# back to right of lens -thickness = 3.0mm -size = 32.5mm 5.0mm -material = "aluminum" -location = "absolute" -offset = 18.75mm 0 -3.5mm - -[support] -# back to left of lens -thickness = 3.0mm -size = 32.5mm 5.0mm -material = "aluminum" -location = "absolute" -offset = -18.75mm 0 -3.5mm - -[support] -# siding, right -thickness = 4.0mm -size = 4.0mm 150mm -material = "aluminum" -location = "absolute" -offset = 33.0mm -50mm 0 - -[support] -# siding, left -thickness = 4.0mm -size = 4.0mm 150mm -material = "aluminum" -location = "absolute" -offset = -33.0mm -50mm 0 - -[support] -# siding, top -thickness = 4.0mm -size = 62mm 4mm -material = "aluminum" -location = "absolute" -offset = 0 23mm 0 - -[support] -# siding, bottom -thickness = 4.0mm -size = 62mm 4mm -material = "aluminum" -location = "absolute" -offset = 0 -123mm 0 diff --git a/htc_wildfire/htc_wildfire_shielded_replace.conf b/htc_wildfire/htc_wildfire_shielded_replace.conf index 4e10149..6cf4f18 100644 --- a/htc_wildfire/htc_wildfire_shielded_replace.conf +++ b/htc_wildfire/htc_wildfire_shielded_replace.conf @@ -3,6 +3,7 @@ type = "monolithic" number_of_pixels = 2592 1944 pixel_size = 0.9um 0.9um sensor_thickness = {}um +#sensor_excess_direction = 'bottom' #sensor_excess = 0.5mm chip_thickness = 1mm @@ -13,84 +14,3 @@ chip_thickness = 1mm # material = "lead" # location = "absolute" # offset = 0 0 -20mm - -# Build the phone layout -[support] -# front screen -thickness = 3.0mm -size = 70mm 150mm -material = "plexiglass" -location = "absolute" -offset = 0 -50mm +3.5mm - -[support] -# lens -thickness = 4.0mm -size = 5mm 5mm -material = "plexiglass" -location = "absolute" -offset = 0 0 -3.5mm - -[support] -# back below lens -thickness = 3.0mm -size = 70mm 122.5mm -material = "aluminum" -location = "absolute" -offset = 0 -63.75mm -3.5mm - -[support] -# back above lens -thickness = 3.0mm -size = 70mm 22.5mm -material = "aluminum" -location = "absolute" -offset = 0 13.75mm -3.5mm - -[support] -# back to right of lens -thickness = 3.0mm -size = 32.5mm 5.0mm -material = "aluminum" -location = "absolute" -offset = 18.75mm 0 -3.5mm - -[support] -# back to left of lens -thickness = 3.0mm -size = 32.5mm 5.0mm -material = "aluminum" -location = "absolute" -offset = -18.75mm 0 -3.5mm - -[support] -# siding, right -thickness = 4.0mm -size = 4.0mm 150mm -material = "aluminum" -location = "absolute" -offset = 33.0mm -50mm 0 - -[support] -# siding, left -thickness = 4.0mm -size = 4.0mm 150mm -material = "aluminum" -location = "absolute" -offset = -33.0mm -50mm 0 - -[support] -# siding, top -thickness = 4.0mm -size = 62mm 4mm -material = "aluminum" -location = "absolute" -offset = 0 23mm 0 - -[support] -# siding, bottom -thickness = 4.0mm -size = 62mm 4mm -material = "aluminum" -location = "absolute" -offset = 0 -123mm 0 diff --git a/htc_wildfire/source_measurement.conf b/htc_wildfire/source_measurement.conf index 09aafe2..30b381e 100644 --- a/htc_wildfire/source_measurement.conf +++ b/htc_wildfire/source_measurement.conf @@ -1,24 +1,24 @@ [AllPix] log_level = "WARNING" log_format = "DEFAULT" -number_of_events = 10 +number_of_events = 100 detectors_file = "detector.conf" model_paths = ./ # Set surrounding material to air [GeometryBuilderGeant4] -world_material = "air" +world_material = "vacuum" # Beam parameters and typical GEANT4 physics lists [DepositionGeant4] physics_list = FTFP_BERT_LIV #enable_pai = true charge_creation_energy = 3.62eV -particle_type = "mu+" -source_energy = 1GeV -source_position = 0 0 -20mm +#particle_type = "mu+" +#source_energy = 1GeV +source_position = 0 0 13.15um source_type = "beam" -beam_direction = 0 0 1 +beam_direction = 1.2246467991473532e-16 0.0 -1.0 beam_divergence = 0.1mrad 0.1mrad max_step_length = 1um output_plots = false @@ -33,7 +33,7 @@ bias_voltage=30V [GenericPropagation] temperature = 293K propagate_electrons = true -propagate_holes = false +propagate_holes = true charge_per_step = 1 integration_time = 10ms output_plots = false @@ -48,11 +48,11 @@ max_depth_distance = 1.0um output_plots = false threshold = 1e threshold_smearing = 1e -adc_smearing = 1e -adc_slope = 25e +qdc_smearing = 1e +qdc_slope = 25e gain = 1.0 electronics_noise = 1e -allow_zero_adc = true +allow_zero_qdc = true # What to write to output # [ROOTObjectWriter] @@ -61,7 +61,6 @@ allow_zero_adc = true [TextWriter] include = "PixelHit" -file_name = "depletion_systematics/100um_mu+_10GeV.txt" # Visualization stuff #[VisualizationGeant4] diff --git a/htc_wildfire/source_measurement_replace.conf b/htc_wildfire/source_measurement_replace.conf index e9a312c..b302086 100644 --- a/htc_wildfire/source_measurement_replace.conf +++ b/htc_wildfire/source_measurement_replace.conf @@ -7,18 +7,18 @@ model_paths = ./ # Set surrounding material to air [GeometryBuilderGeant4] -world_material = "air" +world_material = "vacuum" # Beam parameters and typical GEANT4 physics lists [DepositionGeant4] physics_list = FTFP_BERT_LIV #enable_pai = true charge_creation_energy = 3.62eV -particle_type = "mu+" -source_energy = 1GeV -source_position = 0 0 -20mm +#particle_type = "mu+" +#source_energy = 1GeV +source_position = {} source_type = "beam" -beam_direction = 0 0 1 +beam_direction = {} beam_divergence = 0.1mrad 0.1mrad max_step_length = 1um output_plots = false @@ -33,10 +33,10 @@ bias_voltage=30V [GenericPropagation] temperature = 293K propagate_electrons = true -propagate_holes = false +propagate_holes = true charge_per_step = 1 integration_time = 10ms -output_plots = false +output_plots = {} # combines single charges to sets of charges for processing [SimpleTransfer] @@ -48,11 +48,11 @@ max_depth_distance = 1.0um output_plots = false threshold = 1e threshold_smearing = 1e -adc_smearing = 1e -adc_slope = 25e +qdc_smearing = 1e +qdc_slope = 25e gain = 1.0 electronics_noise = 1e -allow_zero_adc = true +allow_zero_qdc = true # What to write to output # [ROOTObjectWriter] @@ -61,7 +61,6 @@ allow_zero_adc = true [TextWriter] include = "PixelHit" -file_name = "depletion_systematics/100um_mu+_10GeV.txt" # Visualization stuff #[VisualizationGeant4] @@ -70,3 +69,8 @@ file_name = "depletion_systematics/100um_mu+_10GeV.txt" ###simple_view = false #transparency = 0.3 #trajectories_color_mode = "charge" + +# The object writer listens to all output data +#[ROOTObjectWriter] +# specify the output file (default file name is used if omitted) +#file_name = "data.root"