diff --git a/BurstSensitivities.py b/BurstSensitivities.py new file mode 100644 index 0000000..8edf937 --- /dev/null +++ b/BurstSensitivities.py @@ -0,0 +1,221 @@ +#!/usr/bin/env python # +# # +# Autor: Henrike Fleischhack, CUA/NASA-GSFC. # +# # +#------------------------------------------------------------------------------# + +""" +Calculate integral flux sensitivities for bursts at a given time scale and energy range, given effective area and background rates. No ARM cuts! + +This is similar to running MEGALib's SensitivityOptimizer with very open cuts + (and no actual optimization), however: + * Here, we get the background rates from (binned) histograms rather than counting events, + so the resolution is limited by the width of the energy bins. + * Here, we find the effective area by interpolating the effective area vs true energy + curves from mono-energetic simulations and finding the weighted average + given a certain energy spectrum. + The SensitivityOptimizer counts events falling in a given *reconstructed* energy range. + +Note that no cuts are applied here, including no ARM cuts. +This is appropriate for very short time scales (bursts, giant magnetar flares) +but not for the survey! + +usage: BurstSensitivities.py [-h] -f [EAFILES [EAFILES ...]] -b BGFILE + [-em EMIN] [-eM EMAX] [-t EXPOSURE] + [-s SIGNIFICANCE] [-x INDEX] [-c CUTOFF] [-o] + +optional arguments: + -h, --help show this help message and exit + -f [EAFILES [EAFILES ...]], --eafiles [EAFILES [EAFILES ...]] + input .txt file(s) with effective areas (from + plotFigureOfMwerit.py) (default: None) + -b BGFILE, --bgfile BGFILE + input .root file(s) with background rate histograms + (already normalized by exposure time) (default: None) + -em EMIN, --emin EMIN + Minimum energy (MeV) (default: 0.1) + -eM EMAX, --emax EMAX + Maximum energy (MeV) (default: 10) + -t EXPOSURE, --exposure EXPOSURE + Desired exposure time (s) (default: 1) + -s SIGNIFICANCE, --significance SIGNIFICANCE + Desired signal/noise threshold (sigma) (default: 6) + -x INDEX, --index INDEX + Spectral index (default: -2) + -c CUTOFF, --cutoff CUTOFF + Cutoff energy (MeV) (default: inf) + +""" + +import os +import re +import argparse +import numpy as np +import matplotlib.pyplot as plt +import ROOT +from scipy.interpolate import interp1d +from scipy.optimize import minimize, root_scalar + +__description__ = 'Integral flux sensitivities for bursts' + + +"""Command-line switches. +""" +formatter = argparse.ArgumentDefaultsHelpFormatter +PARSER = argparse.ArgumentParser(description=__description__, + formatter_class=formatter) +PARSER.add_argument('-f', '--eafiles', type=str, required=True, nargs='*', + help='input .txt file(s) with effective areas (from plotFigureOfMwerit.py)') +PARSER.add_argument('-b', '--bgfile', type=str, required=True, + help='input .root file(s) with background rate histograms (already normalized by exposure time)') +PARSER.add_argument('-em', '--emin', type=float, default=0.1, + help='Minimum energy (MeV)') +PARSER.add_argument('-eM', '--emax', type=float, default=10, + help='Maximum energy (MeV)') +PARSER.add_argument('-t', '--exposure', type=float, default=1, help='Desired exposure time (s)') +PARSER.add_argument('-s', '--significance', type=float, default=6, help='Desired signal/noise threshold (sigma)') +PARSER.add_argument('-x', '--index', type=float, default=-2, help='Spectral index') +PARSER.add_argument('-c', '--cutoff', type=float, default=np.inf, help='Cutoff energy (MeV)') + +FLAG_STR = '(?<=\_)\w+(?=\.txt)' + + +def get_bg_counts(root_file, emin, emax, type): + """read in ROOT file with BG rate histograms (dN/dE/dT) and return + the rate dN/dT for a given energy range and event type (UC, TC, P)""" + + bgRate = {} + _file0 = ROOT.TFile(root_file) + + + try: + + totalBG = _file0.Get(f"totalBG_{type}") + bin1 = totalBG.FindBin( emin*1e3 ) #energies are in keV! + bin2 = totalBG.FindBin( emax*1e3 ) + rate = totalBG.Integral(bin1, bin2, "width" ) + + except: + rate = None + + return rate + + +def read_effective_areas(ea_files): + """read in txt files with effective area vs true energy (as in compareFigureOfMerit.py) + return a dictionary with the effective areas (as interp1d objects) and event type (UC, TC, P)""" + + theEAs = {} + + for ea_file in ea_files: + + print('Parsing %s ...'%ea_file) + + file_basename = os.path.basename(ea_file) + m = re.search(r'%s'%FLAG_STR, file_basename) + type = m.group(0) + + en, ea = np.loadtxt(ea_file).T + + ea[np.isnan(ea)] = 0 + + f = interp1d(en, ea, kind='cubic') + + theEAs[type] = f + + return theEAs + + +def average_effective_area( EA, emin, emax, index, cutoff): + """return effective area averaged over a given energy range, weighted by a given energy spectrum (cutoff powerlaw)""" + + E = np.linspace( emin, emax, 1000) + dNdE = E**args.index * np.exp(-E/args.cutoff) + + x1 = 0 + x2 = 0 + + for e, n in zip(E, dNdE): + x1 += n * EA(e) + x2 += n + + if x2 == 0: + return 0 + else: + return x1/x2 + + +def get_sensitivity(time, significance, Aeff, BGrate): + """Calculate minimum detectable flux (in photons/cm2/s) for a given exposure time, detection threshold, + effective area, and background rate""" + + arg = np.sqrt((significance**4/4.)+(significance**2*BGrate*time)) + num = (significance**2/2+arg)/(Aeff*time) + return num + +def get_significance(flux, time, Aeff, BGrate): + """Calculate expected significance for a given flux (ph/cm2/s), exposure time, + effective area, and background rate""" + + S = flux * time * Aeff + B = time * BGrate + return S / np.sqrt(S+B) + + + + +if __name__ == '__main__': + + args = PARSER.parse_args() + + theEAs = read_effective_areas(args.eafiles) + + print(f"{args.significance}σ sensitvities for T = {args.exposure:.3g}s") + + Emin = args.emin + Emax = args.emax + + EAs_for_combo = {} + BGs_for_combo = {} + + for type in theEAs.keys(): + + Aeff = average_effective_area( theEAs[type], Emin, Emax, args.index, args.cutoff ) + BGrate = get_bg_counts(args.bgfile, Emin, Emax, type ) + + #only try to calculate the sensitivity if we have events + if BGrate is not None and (Aeff > 0): + sensi = get_sensitivity(args.exposure, args.significance, Aeff, BGrate) + + print( f"{type:2s} Events: {Emin:5.2f} - {Emax:5.2f} MeV, " + f" Aeff = {Aeff:6.4g} cm2, " + f"BG = {BGrate:5.1f} Hz, min. Flux = {sensi:7.4g} ph/cm2/s" ) + + EAs_for_combo[type] = Aeff + BGs_for_combo[type] = BGrate + + def get_combined_significance(flux): + + sig = 0 + for type in EAs_for_combo.keys(): + sig += get_significance( flux, args.exposure, EAs_for_combo[type], BGs_for_combo[type])**2 + + return np.sqrt(sig) + + def function_to_root( flux ): + return get_combined_significance(flux) - args.significance + + solution = root_scalar(function_to_root, x0 = 0, x1 = sensi) + min_flux = solution.root + + out_string = f"Combo: min. Flux = {min_flux:7.4g} ph/cm2/s, significance = " + total_sig = 0 + for type in EAs_for_combo.keys(): + sig = get_significance( min_flux, args.exposure, EAs_for_combo[type], BGs_for_combo[type]) + out_string += f" {sig:.3g}[{type}] ⊕" + total_sig+=sig**2 + out_string = out_string[:-2] + total_sig = np.sqrt(total_sig) + out_string += f" = {total_sig:.3g}" + + print( out_string ) diff --git a/EventAnalysis.py b/EventAnalysis.py index c25dfd3..8870fe8 100644 --- a/EventAnalysis.py +++ b/EventAnalysis.py @@ -55,6 +55,7 @@ import scipy.optimize import math import glob +import gzip try: import matplotlib.pyplot as plot @@ -339,6 +340,10 @@ def parse(filename, sourceTheta=1.0, testnum=-1): qualityOfComptonReconstruction = [] qualityOfPairReconstruction = [] + # Original (true) information + true_direction_PairEvents = [] + true_energy_PairEvents = [] + # Read the number of lines in the file command = 'wc %s' % filename @@ -370,7 +375,13 @@ def parse(filename, sourceTheta=1.0, testnum=-1): # Loop through the .tra file print('\nParsing: %s' % filename) - for line in fileinput.input([filename]): + + if filename[-2:] == "gz": #compressed file + lines = gzip.open(filename, mode="rt") + else: + lines = fileinput.input([filename]) + + for line in lines: try: sys.stdout.write("Progress: %d%% \r" % (lineNumber/totalNumberOfLines * 100) ) @@ -390,6 +401,10 @@ def parse(filename, sourceTheta=1.0, testnum=-1): if 'ET ' in line: + # new event -- reset true information just to be safe + trueDirection = None + trueEnergy = None + # Split the line lineContents = line.split() @@ -429,6 +444,15 @@ def parse(filename, sourceTheta=1.0, testnum=-1): # In pair, it's mostly 'TooManyHistInCSR' numberOfBadEvents = numberOfBadEvents + 1 + ####### True Information ####### + + if 'OI' in line and skipEvent == False: + + lineContents = line.split() + xTrue, yTrue, zTrue = lineContents[4:7] + trueDirection = [-float(xTrue), -float(yTrue), -float(zTrue)] + trueEnergy = float( lineContents[10] ) + ####### Compton Events ####### # Extract the Compton event energy information @@ -547,8 +571,14 @@ def parse(filename, sourceTheta=1.0, testnum=-1): # Calculate the vector between the second and first interaction directionVector2 = position2 - position1 - # Calculate the vector between the first interaction and the origin of the original gamma-ray - directionVector1 = position1 - position0 + if trueDirection is None: + + # Calculate the vector between the first interaction and the origin of the original gamma-ray + directionVector1 = position1 - position0 + + else: + directionVector1 = [-x for x in trueDirection] + # Calculate the product of the vector magnitudes product = numpy.linalg.norm(directionVector1) * numpy.linalg.norm(directionVector2) @@ -599,6 +629,10 @@ def parse(filename, sourceTheta=1.0, testnum=-1): # Save the position position_pairConversion.append([x1,y1,z1]) + # True information should have been read in earlier + true_direction_PairEvents.append( trueDirection ) + true_energy_PairEvents.append( trueEnergy ) + # Extract the pair electron information @@ -607,36 +641,38 @@ def parse(filename, sourceTheta=1.0, testnum=-1): # Split the line lineContents = line.split() + if len(lineContents) >= 6: #PE token can be used for either pair events or photoelectric effect. + + # Get the electron information + energy_pairElectron.append(float(lineContents[1])) + energy_pairElectron_error.append(float(lineContents[2])) - # Get the electron information - energy_pairElectron.append(float(lineContents[1])) - energy_pairElectron_error.append(float(lineContents[2])) - - # Get the direction of the pair electron - x = float(lineContents[3]) - y = float(lineContents[4]) - z = float(lineContents[5]) + # Get the direction of the pair electron + x = float(lineContents[3]) + y = float(lineContents[4]) + z = float(lineContents[5]) - # Store the direction of the pair electron - direction_pairElectron.append([x,y,z]) + # Store the direction of the pair electron + direction_pairElectron.append([x,y,z]) # Extract the pair positron information if 'PP ' in line and skipEvent == False: # Split the line lineContents = line.split() + if len(lineContents) >= 6: #PP token can be used for either pair events or photoelectric effect. - # Get the electron information - energy_pairPositron.append(float(lineContents[1])) - energy_pairPositron_error.append(float(lineContents[2])) + # Get the electron information + energy_pairPositron.append(float(lineContents[1])) + energy_pairPositron_error.append(float(lineContents[2])) - # Get the direction of the pair electron - x = float(lineContents[3]) - y = float(lineContents[4]) - z = float(lineContents[5]) + # Get the direction of the pair electron + x = float(lineContents[3]) + y = float(lineContents[4]) + z = float(lineContents[5]) - # Store the direction of the pair electron - direction_pairPositron.append([x,y,z]) + # Store the direction of the pair electron + direction_pairPositron.append([x,y,z]) if 'PI ' in line and skipEvent == False: @@ -710,24 +746,29 @@ def parse(filename, sourceTheta=1.0, testnum=-1): events['qualityOfPairReconstruction'] = numpy.array(qualityOfPairReconstruction).astype(float) events['time'] = numpy.array(time).astype(float) events['deltime'] = numpy.append(0.,events['time'][1:]-events['time'][0:len(events['time'])-1]) + events['true_direction_PairEvents'] = true_direction_PairEvents + events['true_energy_PairEvents'] = true_energy_PairEvents + events['numberOfPhotoElectricEffectEvents'] = numberOfPhotoElectricEffectEvents # Print some event statistics print("\n\nStatistics of Event Selection") print("***********************************") print("Total number of analyzed events: %s" % (numberOfComptonEvents + numberOfPairEvents)) - if numberOfComptonEvents + numberOfPairEvents == 0: + if numberOfComptonEvents + numberOfPairEvents + numberOfPhotoElectricEffectEvents == 0: print("No events pass selection") events=False return events + numberOfTotalEvents = numberOfComptonEvents + numberOfPairEvents + numberOfUnknownEventTypes + numberOfPhotoElectricEffectEvents print("") - print("Number of unknown events: %s (%i%%)" % (numberOfUnknownEventTypes, 100*numberOfUnknownEventTypes/(numberOfComptonEvents + numberOfPairEvents + numberOfUnknownEventTypes))) - print("Number of pair events: %s (%i%%)" % (numberOfPairEvents, 100*numberOfPairEvents/(numberOfComptonEvents + numberOfPairEvents + numberOfUnknownEventTypes))) - print("Number of Compton events: %s (%i%%)" % (numberOfComptonEvents, 100*numberOfComptonEvents/(numberOfComptonEvents + numberOfPairEvents + numberOfUnknownEventTypes))) + print("Number of unknown events: %s (%i%%)" % (numberOfUnknownEventTypes, 100*numberOfUnknownEventTypes/numberOfTotalEvents)) + print("Number of pair events: %s (%i%%)" % (numberOfPairEvents, 100*numberOfPairEvents/numberOfTotalEvents)) + print("Number of Compton events: %s (%i%%)" % (numberOfComptonEvents, 100*numberOfComptonEvents/numberOfTotalEvents)) if numberOfComptonEvents > 0: print(" - Number of tracked electron events: %s (%i%%)" % (numberOfTrackedElectronEvents, 100.0*(float(numberOfTrackedElectronEvents)/numberOfComptonEvents))) print(" - Number of untracked electron events: %s (%i%%)" % (numberOfUntrackedElectronEvents, 100*(float(numberOfUntrackedElectronEvents)/numberOfComptonEvents))) + print("Number of p.e. events: %s (%i%%)" % (numberOfPhotoElectricEffectEvents, 100*numberOfPhotoElectricEffectEvents/numberOfTotalEvents)) print("") print("") @@ -795,14 +836,14 @@ def getARMForComptonEvents(events, numberOfBins=100, phiRadius=10, onlyTrackedEl exit() # Calculate the Compton scattering angle - value = 1 - electron_mc2 * (1/energy_firstScatteredPhoton - 1/(energy_recoiledElectron + energy_firstScatteredPhoton)); + value = 1 - electron_mc2 * (1/energy_firstScatteredPhoton - 1/(energy_recoiledElectron + energy_firstScatteredPhoton)) # Keep only sane results # index = numpy.where( (value > -1) | (value < 1) ) # value = value[index] # Calculate the theoretical phi angle (in radians) - phi_Theoretical = numpy.arccos(value); + phi_Theoretical = numpy.arccos(value) # Calculate the difference between the tracker reconstructed scattering angle and the theoretical scattering angle dphi = numpy.degrees(phi_Tracker) - numpy.degrees(phi_Theoretical) @@ -817,12 +858,12 @@ def getARMForComptonEvents(events, numberOfBins=100, phiRadius=10, onlyTrackedEl gs = gridspec.GridSpec(4,1) ax1 = plot.subplot(gs[:3, :]) - if len(dphi[selection]) < 10: - print("The number of events is not sufficient (<10)") + if len(dphi[selection]) < 50: + print("The number of events is not sufficient (<50)") return numpy.nan, numpy.nan - elif len(dphi[selection]) < 500: - print("The number of events is not sufficient, so that 'numberOfBins' and 'phiRadius' changed.") - numberOfBins = 25 + #elif len(dphi[selection]) < 500: + # print("The number of events is not sufficient, so that 'numberOfBins' and 'phiRadius' changed.") + # numberOfBins = 25 # Create the histogram histogram_angleResults = ax1.hist(dphi[selection], numberOfBins, color='#3e4d8b', alpha=0.9, histtype='stepfilled') @@ -830,11 +871,11 @@ def getARMForComptonEvents(events, numberOfBins=100, phiRadius=10, onlyTrackedEl # Extract the binned data and bin locations dphi_binned = histogram_angleResults[0] - if max(dphi_binned) < 10: - numberOfBins = 20 - histogram_angleResults = ax1.hist(dphi[selection], numberOfBins, color='#3e4d8b', alpha=0.9, histtype='stepfilled') - ax1.set_xlim([-1*phiRadius,phiRadius]) - dphi_binned = histogram_angleResults[0] + #if max(dphi_binned) < 10: + # numberOfBins = 20 + # histogram_angleResults = ax1.hist(dphi[selection], numberOfBins, color='#3e4d8b', alpha=0.9, histtype='stepfilled') + # ax1.set_xlim([-1*phiRadius,phiRadius]) + # dphi_binned = histogram_angleResults[0] bins = histogram_angleResults[1] bincenters = 0.5*(bins[1:]+bins[:-1]) @@ -1022,14 +1063,23 @@ def getARMForPairEvents(events, sourceTheta=0, numberOfBins=100, angleFitRange=[ # Get the position of the gamma conversion position_conversion = events['position_pairConversion'][index] + direction_source = events['true_direction_PairEvents'][index] + #print( direction_source ) + + if direction_source is None: - # Get the x-axis offset based on the theta of the source. This assumes phi=0 - # Note: For Compton events adjusting the theta of the source happens in the parser - # for pair events, it happens here in the ARM calculation. - dx = numpy.tan(numpy.radians(sourceTheta)) * (position_conversion[2] - dz) + # Try to reconstruct the source direction from the cosTheta in the file name. + # Get the x-axis offset based on the theta of the source. This assumes phi=0 + # Note: For Compton events adjusting the theta of the source happens in the parser + # for pair events, it happens here in the ARM calculation. + dx = numpy.tan(numpy.radians(sourceTheta)) * (position_conversion[2] - dz) - # Set the origin position of the original gamma-ray - position_source = [position_conversion[0]-dx, position_conversion[1], dz] + # Set the origin position of the original gamma-ray + position_source = [position_conversion[0]-dx, position_conversion[1], dz] + + # Calculate the vector between the first interaction and the origin of the original gamma-ray + direction_source = -1*(position_conversion - position_source) + #print( direction_source ) # Get the electron and positron direction vectors. These are unit vectors. direction_electron = events['direction_pairElectron'][index] @@ -1046,9 +1096,6 @@ def getARMForPairEvents(events, sourceTheta=0, numberOfBins=100, angleFitRange=[ # Invert the bisect vector to obtain the reconstructed source vector direction_source_reconstructed = -1*direction_bisect - # Calculate the vector between the first interaction and the origin of the original gamma-ray - direction_source = -1*(position_conversion - position_source) - # Calculate the distance of the conversion point to the top-center of the spacecraft position_topCenter = numpy.array([0,0,60]) distance = numpy.linalg.norm(position_conversion-position_topCenter) @@ -1420,7 +1467,7 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N if len(energy_pairReconstructedPhoton[selection]) < 10: print("The number of events is too small (<10)") - return numpy.nan, numpy.nan, numpy.nan + return numpy.nan, numpy.nan, numpy.nan, numpy.nan, numpy.nan, numpy.nan elif len(energy_pairReconstructedPhoton[selection]) < 100: numberOfBins = 30 @@ -1455,7 +1502,7 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N y_fit = skewedGaussian(bincenters, optimizedParameters[0], optimizedParameters[1], optimizedParameters[2], optimizedParameters[3]) except Exception as msg: print("{}".format(msg)) - return numpy.nan, numpy.nan, numpy.nan + return numpy.nan, numpy.nan, numpy.nan, numpy.nan, numpy.nan, numpy.nan # Get the max of the fit fitMax = bincenters[numpy.argmax(y_fit)] @@ -1466,18 +1513,8 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N FWHM = x2-x1 # Approximate sigma as FWHM/2 - sigma=FWHM/2. + sigma=FWHM/2.355 - # Print some statistics - print("\n\nStatistics of energy histogram and fit (pair events)") - print("********************************************************") - print("Number of Compton and pair events in histogram: %s (%s%%)" % ( len(energy_pairReconstructedPhoton[selection]), 100*len(energy_pairReconstructedPhoton[selection])/(len(energy_pairReconstructedPhoton)) )) - print("") - print("Fitting in range: ", energyFitRange[0], energyFitRange[1]) - print("Max of fit: %s keV" % fitMax) - print("FWHM of fit: %s keV" % FWHM) - print("sigma of fit: %s keV" % sigma) - print("") # Set the plot size @@ -1499,10 +1536,7 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N ax1.plot(bincenters, y_fit, color='darkred', linewidth=2) ax1.plot([x1,x2],[numpy.max(y_fit)/2.,numpy.max(y_fit)/2.], color='darkred', linestyle='--', linewidth=2) - # Annotate the plot - ax1.text(0.03, 0.9, "Max = %.3f keV\nFWHM = %.3f keV" % (fitMax, FWHM), verticalalignment='bottom', horizontalalignment='left', transform=ax1.transAxes, color='black', fontsize=12) - - # Create a subplot for the residuals + # Create a subplot for the residuals ax2 = plot.subplot(gs[3, :]) # Plot the residuals @@ -1515,6 +1549,19 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N ax1.xaxis.set_minor_locator(AutoMinorLocator(4)) ax1.yaxis.set_minor_locator(AutoMinorLocator(4)) ax2.xaxis.set_minor_locator(AutoMinorLocator(4)) + + selection2 = numpy.where( qualityOfPairReconstruction <= qualityCut ) + low, med, up = numpy.quantile( energy_pairReconstructedPhoton[selection2], ( .16, 0.5, .84) ) + + ax1.axvline( low, color="0.5", linestyle=":", linewidth=1.5) + ax1.axvline( med, color="0.5", linestyle=":", linewidth=3) + ax1.axvline( up, color="0.5", linestyle=":", linewidth=1.5) + + # Annotate the plot + ax1.text(0.03, 0.9, "Max = %.3f keV\nFWHM = %.3f keV\nmedian = %.3f\n68%% containment half width = %3f keV" % (fitMax, FWHM, med, (up-low)), verticalalignment='bottom', horizontalalignment='left', transform=ax1.transAxes, color='black', fontsize=12) + + + plot.savefig(fileBase+"_energyResolution_pairs.png") # Show the plot if showPlots == True: @@ -1522,7 +1569,21 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N else: plot.close() - return fitMax, FWHM, sigma + # Print some statistics + print("\n\nStatistics of energy histogram and fit (pair events)") + print("********************************************************") + print("Number of Compton and pair events in histogram: %s (%s%%)" % ( len(energy_pairReconstructedPhoton[selection]), 100*len(energy_pairReconstructedPhoton[selection])/(len(energy_pairReconstructedPhoton)) )) + print("") + print("Fitting in range: ", energyFitRange[0], energyFitRange[1]) + print("Median reco energy: %s keV" % med) + print("68%% Containment Interval: %s - %s keV" % (low, up) ) + print("Max of fit: %s keV" % fitMax) + print("FWHM of fit: %s keV" % FWHM) + print("sigma of fit: %s keV" % sigma) + print("") + + + return fitMax, FWHM, sigma, low, med, up ########################################################################################## @@ -2022,7 +2083,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles # Check to see if the user supplied a directory. If so, include all .tra files in the directory if directory != None: - filenames = glob.glob(directory + '/*.tra') + filenames = glob.glob( '{}/*.tra'.format(directory)) + glob.glob( '{}/*.tra.gz'.format(directory)) # Check if the user supplied a single file vs a list of files @@ -2067,18 +2128,30 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles # Loop through the user specified filename(s) and extract the energy and angular resolution measurements for filename, energy, angle in zip(filenames, energies, angles): + output_filename = filename.replace('.tra','.log').replace(".gz", "") + if os.path.exists(output_filename): + print("Skipping: %s %s Cos %s %s" % (energy, energySearchUnit, angle, filename)) + continue + if parsing: print("Parsing: %s %s Cos %s %s" % (energy, energySearchUnit, angle, filename)) # Parse the .tra file obtained from revan events = parse(filename, sourceTheta=angle) + if not events: #no events pass the selection + continue # Calculate the source theta in degrees source_theta = numpy.arccos(angle)*180./numpy.pi + # Don't bother measuring the energy and angular resolutuon values for Compton events above the specified maximumComptonEnergy if energy <= maximumComptonEnergy: - if energy >= 3: - phiRadiusCompton = phiRadiusCompton/3. + + phiRadiusCompton = 15 + if energy < 0.3: + phiRadiusCompton = 30 + if energy < 0.2: + phiRadiusCompton = 60 print("--------- All Compton Events ---------") # Calculate the energy resolution for Compton events @@ -2146,7 +2219,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles print("\n\nCalculating the energy resolution for pair events...") print("EventAnalysis.getEnergyResolutionForPairEvents(events, numberOfBins=100)") fileBase = "%sMeV_Cos%s" % (energy,angle) - fitMax, FWHM_pairComptonEvents, sigma_pair = getEnergyResolutionForPairEvents(events, numberOfBins=100, energyFitRange=True, showPlots=showPlots,fileBase=fileBase) + fitMax_pair, FWHM_pairComptonEvents, sigma_pair, low_pair, med_pair, up_pair = getEnergyResolutionForPairEvents(events, numberOfBins=100, energyFitRange=True, showPlots=showPlots,fileBase=fileBase) # Calculate the angular resolution measurement (ARM) for pair events print("\n\nCalculating the angular resolution measurement for pair events...") @@ -2154,12 +2227,15 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles angles, openingAngles, contaimentData_68, contaimentBinned_68 = getARMForPairEvents(events, openingAngleMax=openingAngleMax, sourceTheta=source_theta, numberOfBins=100, showDiagnosticPlots=False, showPlots=showPlots, filename=filename) else: + fitMax_pair = numpy.nan + med_pair = numpy.nan + up_pair = numpy.nan + low_pair = numpy.nan sigma_pair = numpy.nan contaimentData_68 = numpy.nan FWHM_pairComptonEvents = numpy.nan # Open the results filename for writing - output_filename = filename.replace('.tra','.log') output = open(output_filename, 'w') # Write the results to disk @@ -2182,9 +2258,13 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles output.write("Pair Events Reconstructed: %s\n" % events['numberOfPairEvents']) output.write("Pair Energy Resolution (keV): %s\n" % sigma_pair) #FWHM_pairComptonEvents + output.write("Pair Energy FitMax (keV): %s\n" % fitMax_pair) + output.write("Pair Energy Median (keV): %s\n" % med_pair) + output.write("Pair Energy 68%% Containment Interval (keV): %s\n" % (up_pair - low_pair) ) output.write("Pair Angular Containment (68%%): %s\n" % contaimentData_68) output.write("Events Not Reconstructed Flagged as Bad: %s\n" % events['numberOfBadEvents']) + output.write("Number of Photo Electric Effect Events: %s\n" % events['numberOfPhotoElectricEffectEvents']) # Close the file output.close() @@ -2216,7 +2296,7 @@ def getTriggerEfficiency(filename=None, directory=None, save=True, savefile=None # Check to see if the user supplied a directory. If so, include all .sim files in the directory if directory != None: print("\nSearching: %s\n" % directory) - filenames = glob.glob(directory + '/*.sim') + filenames = glob.glob(directory + '/*.sim' ) + glob.glob(directory + '/*.sim.gz' ) # Check if the user supplied a single file vs a list of files if isinstance(filename, list) == False and filename != None: @@ -2241,7 +2321,10 @@ def getTriggerEfficiency(filename=None, directory=None, save=True, savefile=None lookback = 1000 IDs = [] while len(IDs) == 0: - command = "tail -n %d %s" % (lookback, filename) + if filename[-2:] == "gz": + command = "zcat %s | tail -n %d" % (filename, lookback) + else: + command = "tail -n %d %s" % (lookback, filename) output = os.popen(command).read() IDs = [line for line in output.split('\n') if "ID" in line] lookback += 1000 @@ -2319,7 +2402,7 @@ def getRevanTriggerEfficiency(filename=None, directory=None, save=True, savefile # Check to see if the user supplied a directory. If so, include all .tra files in the directory if directory != None: print("\nSearching: %s\n" % directory) - filenames = glob.glob(directory + '/*.tra') + filenames = glob.glob( './*.tra') + glob.glob( './*.tra.gz') # Check if the user supplied a single file vs a list of files if isinstance(filename, list) == False and filename != None: diff --git a/FigureOfMeritPlotter.py b/FigureOfMeritPlotter.py index 660e421..dbc477d 100644 --- a/FigureOfMeritPlotter.py +++ b/FigureOfMeritPlotter.py @@ -258,18 +258,12 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal # Get the filename simulationName = lineContents[0].split('/')[-1] - # Create a key for this simulation name - data[simulationName] = [] - # Get the number of triggered and simulated events numberOfTriggeredEvents = lineContents[1] numberOfSimulatedEvents = lineContents[2] - # Add the number of simulated to the results dictionary - data[simulationName].append(numberOfSimulatedEvents) - # Generate the log filename - analysisLog = directory + '/' + simulationName.replace('.sim', '.log') + analysisLog = directory + '/' + simulationName.replace('.sim', '.log').replace(".gz", "") try: if silent==False: @@ -279,9 +273,18 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal except: print("*** Could not find log file: %s" % analysisLog) + continue + + # Create a key for this simulation name + data[simulationName] = [] + + # Add the number of simulated to the results dictionary + data[simulationName].append(numberOfSimulatedEvents) # Loop through the analysis log file numberOfPairEventsIdeal = None + numberOfPhotoElectricEffectEvents = 0 + for analysisLogLine in analysisLogLines: if "Results for simulation:" in analysisLogLine: @@ -330,6 +333,9 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal if "Pair Events Ideal: " in analysisLogLine: numberOfPairEventsIdeal = analysisLogLine.split()[-1] print("pair events ideal: ", numberOfPairEventsIdeal) + + if "Number of Photo Electric Effect Events" in analysisLogLine: + numberOfPhotoElectricEffectEvents = analysisLogLine.split()[-1] # Add all the values to the results dictionary @@ -339,7 +345,7 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal holder=-999. # This is required for historic reasons if float(energy) <=30.0: - data[simulationName].append([holder, holder, holder, holder, FWHM_energyComptonEvents, holder, holder, FWHM_angleComptonEvents, numberOfComptonEvents]) + data[simulationName].append([numberOfPhotoElectricEffectEvents, holder, holder, holder, FWHM_energyComptonEvents, holder, holder, FWHM_angleComptonEvents, numberOfComptonEvents]) data[simulationName].append([holder, holder, holder, holder, FWHM_energyTrackedComptonEvents, holder, holder, FWHM_angleTrackedComptonEvents, numberOfTrackedComptonEvents]) data[simulationName].append([holder, holder, holder, holder, FWHM_energyUntrackedComptonEvents, holder, holder, FWHM_angleUntrackedComptonEvents, numberOfUntrackedComptonEvents]) else: @@ -393,7 +399,7 @@ def plotAngularResolution(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], xlog=Tr energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = float(half.replace('.inc1.id1.sim','')) + angle = float(half.replace('.inc1.id1.sim','').replace(".gz", "")) if angle == angleSelection: Energy.append(energy) @@ -432,19 +438,19 @@ def plotAngularResolution(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], xlog=Tr # writing txt files results_txt_TC.write("# Energy[MeV] AngRes_TkrCompton[deg]\n") for ii, en in enumerate(Energy[i]): - results_txt_TC.write("%.1f\t%.1f\n"%(en, st[i][ii])) + results_txt_TC.write("%.3f\t%.3f\n"%(en, st[i][ii])) results_txt_TC.close() print('Created %s_AngRes_Cos%s_TC.txt ...!'%(txtOutfileLabel, angleSelection)) results_txt_UC.write("# Energy[MeV] AngRes_UntkrCompton[deg]\n") for ii, en in enumerate(Energy[j]): - results_txt_UC.write("%.1f\t%.1f\n"%(en, sut[j][ii])) + results_txt_UC.write("%.3f\t%.3f\n"%(en, sut[j][ii])) results_txt_UC.close() print('Created %s_AngRes_Cos%s_UC.txt ...!'%(txtOutfileLabel, angleSelection)) results_txt_P.write("# Energy[MeV] AngRes_Pair[deg]\n") for ii, en in enumerate(Energy[k]): - results_txt_P.write("%.1f\t%.1f\n"%(en, sp[k][ii])) + results_txt_P.write("%.3f\t%.3f\n"%(en, sp[k][ii])) results_txt_P.close() print('Created %s_AngRes_Cos%s_P.txt ...!'%(txtOutfileLabel, angleSelection)) @@ -562,7 +568,7 @@ def plotAngularResolutionVsAngle(data, energySelections=None, xlog=False, ylog=F energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = numpy.array(float(half.replace('.inc1.id1.sim',''))) + angle = numpy.array(float(half.replace('.inc1.id1.sim','').replace(".gz", ""))) angle = round(numpy.degrees(numpy.arccos(angle))) if energy == energySelection: @@ -688,7 +694,7 @@ def plotEnergyResolution(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], xlog=Tru energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = numpy.array(float(half.replace('.inc1.id1.sim',''))) + angle = numpy.array(float(half.replace('.inc1.id1.sim','').replace(".gz", ""))) if angle == angleSelection: Energy.append(energy) @@ -728,19 +734,19 @@ def plotEnergyResolution(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], xlog=Tru # writing txt files results_txt_TC.write("# Energy[MeV] EnRes_TkrCompton[FWHM/Energy]\n") for ii, en in enumerate(Energy[i]): - results_txt_TC.write("%.1f\t%.3f\n"%(en, st[ii])) + results_txt_TC.write("%.3f\t%.3f\n"%(en, st[ii])) results_txt_TC.close() print('Created %s_EnRes_Cos%s_TC.txt ...!'%(txtOutfileLabel, angleSelection)) results_txt_UC.write("# Energy[MeV] EnRes_UntkrCompton[FWHM/Energy]\n") for ii, en in enumerate(Energy[j]): - results_txt_UC.write("%.1f\t%.3f\n"%(en, sut[ii])) + results_txt_UC.write("%.3f\t%.3f\n"%(en, sut[ii])) results_txt_UC.close() print('Created %s_EnRes_Cos%s_UC.txt ...!'%(txtOutfileLabel, angleSelection)) results_txt_P.write("# Energy[MeV] EnRes_Pair[FWHM/Energy]\n") for ii, en in enumerate(Energy[k]): - results_txt_P.write("%.1f\t%.3f\n"%(en, sp[ii])) + results_txt_P.write("%.3f\t%.3f\n"%(en, sp[ii])) results_txt_P.close() print('Created %s_EnRes_Cos%s_P.txt ...!'%(txtOutfileLabel, angleSelection)) @@ -861,7 +867,7 @@ def plotEnergyResolutionVsAngle(data, energySelections=None, xlog=False, ylog=Fa energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = numpy.array(float(half.replace('.inc1.id1.sim',''))) + angle = numpy.array(float(half.replace('.inc1.id1.sim','').replace(".gz", ""))) if energy == energySelection: Angle.append(angle) @@ -957,7 +963,7 @@ def plotEnergyResolutionVsAngle(data, energySelections=None, xlog=False, ylog=Fa def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False, xlog=True, ylog=False, save=False, show=True, collapse=False, - SurroundingSphere=150, txtOutfileLabel='xxx'): + SurroundingSphere=150, txtOutfileLabel='xxx', includePE = False): if hasattr(angleSelections, '__iter__') == False: angleSelections = [angleSelections] @@ -977,23 +983,27 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False results_txt_TC = open( '%s_Aeff_Cos%s_TC.txt' % (txtOutfileLabel, angleSelection), 'w') results_txt_UC = open( '%s_Aeff_Cos%s_UC.txt' % (txtOutfileLabel, angleSelection), 'w') results_txt_P = open( '%s_Aeff_Cos%s_P.txt' % (txtOutfileLabel, angleSelection), 'w') + results_txt_S = open( '%s_Aeff_Cos%s_S.txt' % (txtOutfileLabel, angleSelection), 'w') Energy = [] EffectiveArea_Tracked = [] EffectiveArea_Untracked = [] EffectiveArea_Pair = [] + EffectiveArea_SingleSite = [] for key in data.keys(): energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = float(half.replace('.inc1.id1.sim','')) + angle = float(half.replace('.inc1.id1.sim','').replace(".gz", "")) #print energy, angle if angle == angleSelection: #print data numberOfSimulatedEvents = float(data[key][0]) + numberOfPhotoElectricEffectEvents = float(data[key][1][0]) + if ideal: # This removes the event selection on the final Aeff calculation # It does not change anything from the FWHM or the 68% containment @@ -1021,12 +1031,14 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False effectiveArea_tracked = (numberOfReconstructedEvents_tracked/numberOfSimulatedEvents) * math.pi * SurroundingSphere**2 effectiveArea_untracked = (numberOfReconstructedEvents_untracked/numberOfSimulatedEvents) * math.pi * SurroundingSphere**2 effectiveArea_pair = (numberOfReconstructedEvents_pair/numberOfSimulatedEvents) * math.pi * SurroundingSphere**2 + effectiveArea_single = (numberOfPhotoElectricEffectEvents/numberOfSimulatedEvents) * math.pi * SurroundingSphere**2 # Store the results Energy.append(energy) EffectiveArea_Tracked.append(effectiveArea_tracked) EffectiveArea_Untracked.append(effectiveArea_untracked) EffectiveArea_Pair.append(effectiveArea_pair) + EffectiveArea_SingleSite.append(effectiveArea_single) # Convert everything to a numpy array @@ -1034,6 +1046,7 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False EffectiveArea_Tracked = numpy.array(EffectiveArea_Tracked) EffectiveArea_Untracked = numpy.array(EffectiveArea_Untracked) EffectiveArea_Pair = numpy.array(EffectiveArea_Pair) + EffectiveArea_SingleSite = numpy.array(EffectiveArea_SingleSite) # Sort by energy i = [numpy.argsort(Energy)] @@ -1041,6 +1054,7 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False EffectiveArea_Tracked = EffectiveArea_Tracked[i] EffectiveArea_Untracked = EffectiveArea_Untracked[i] EffectiveArea_Pair = EffectiveArea_Pair[i] + EffectiveArea_SingleSite = EffectiveArea_SingleSite[i] #EffectiveArea_UntrackedSiStart=numpy.array([96.0/2884228*70685.8, 14936.0/2132319*70685.8, 21777.0/1744552*70685.8,\ # 17861.0/1568899*70685.8, 13722.0/1567836*70685.8, #numpy.nan, numpy.nan, numpy.nan, numpy.nan, numpy.nan, \ @@ -1056,21 +1070,27 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False # writing txt files results_txt_TC.write("# Energy[MeV] Aeff_TkrCompton[cm2]\n") for ii, en in enumerate(Energy): - results_txt_TC.write("%.1f\t%.1f\n"%(en, EffectiveArea_Tracked[ii])) + results_txt_TC.write("%.3f\t%.3f\n"%(en, EffectiveArea_Tracked[ii])) results_txt_TC.close() print('Created %s_Aeff_Cos%s_TC.txt ...!'%(txtOutfileLabel, angleSelections)) results_txt_UC.write("# Energy[MeV] Aeff_UntkrCompton[cm2]\n") for ii, en in enumerate(Energy): - results_txt_UC.write("%.1f\t%.1f\n"%(en, EffectiveArea_Untracked[ii])) + results_txt_UC.write("%.3f\t%.3f\n"%(en, EffectiveArea_Untracked[ii])) results_txt_UC.close() print('Created %s_Aeff_Cos%s_UC.txt ...!'%(txtOutfileLabel, angleSelection)) - + results_txt_P.write("# Energy[MeV] Aeff_Pair[cm2]\n") for ii, en in enumerate(Energy): - results_txt_P.write("%.1f\t%.1f\n"%(en, EffectiveArea_Pair[ii])) + results_txt_P.write("%.3f\t%.3f\n"%(en, EffectiveArea_Pair[ii])) results_txt_P.close() print('Created %s_Aeff_Cos%s_P.txt ...!'%(txtOutfileLabel, angleSelection)) + + results_txt_S.write("# Energy[MeV] Aeff_SingleSite[cm2]\n") + for ii, en in enumerate(Energy): + results_txt_S.write("%.3f\t%.3f\n"%(en, EffectiveArea_SingleSite[ii])) + results_txt_S.close() + print('Created %s_Aeff_Cos%s_S.txt ...!'%(txtOutfileLabel, angleSelection)) # Plot the data @@ -1087,6 +1107,9 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False plot.scatter(Energy, EffectiveArea_Untracked, color='blue') plot.plot(Energy, EffectiveArea_Untracked, color='blue', alpha=0.5, lw=3, label='Untracked Compton') + plot.scatter(Energy, EffectiveArea_SingleSite, color='orange') + plot.plot(Energy, EffectiveArea_SingleSite, color='orange', alpha=0.5, lw=3, label='Single Site') + #plot.scatter(Energy, EffectiveArea_UntrackedSiStart, color='blue') #plot.plot(Energy, EffectiveArea_UntrackedSiStart, color='blue', linestyle="--", alpha=0.5, lw=3, label='Untracked Compton in Silicon') @@ -1170,7 +1193,7 @@ def tabulateEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=F for key in data.keys(): energy = float(key.split('_')[1].replace('MeV', '')) half = key.split('_')[2].replace('Cos', '') - angle = float(half.replace('.inc1.id1.sim', '')) + angle = float(half.replace('.inc1.id1.sim', '').replace(".gz", "")) #print energy, angle @@ -1277,7 +1300,7 @@ def plotEffectiveAreaVsAngle(data, energySelections=None, ideal=False, xlog=Fals energy = float(key.split('_')[1].replace('MeV','')) #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = numpy.array(float(half.replace('.inc1.id1.sim',''))) + angle = numpy.array(float(half.replace('.inc1.id1.sim','').replace(".gz", ""))) angle = round(numpy.degrees(numpy.arccos(angle))) if energy == energySelection: @@ -1399,7 +1422,7 @@ def resultsToFits(data, outfile='output.fits'): for key in data.keys(): energy = float(key.split('_')[1].replace('MeV','')) half = key.split('_')[2].replace('Cos','') - angle = float(half.replace('.inc1.id1.sim','')) + angle = float(half.replace('.inc1.id1.sim','').replace(".gz", "")) if energy not in Energy: Energy.append(energy) if angle not in Angle: @@ -1421,9 +1444,9 @@ def resultsToFits(data, outfile='output.fits'): ########################################################################################## -def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, ideal=False, doPSF=None, \ +def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3*365.25*24*3600*0.2, ideal=False, doPSF=None, \ xlog=True, ylog=True, save=False, doplot=False, showbackground=False, uniterg=False,doRealBkg=True, \ - SurroundingSphere=150): + SurroundingSphere=150, txtOutfileLabel='xxx', data_nominal = None): #background = numpy.array([0.00346008, 0.00447618, 0.00594937, 0.00812853, 0.0100297, 0.0124697, 0.0161290]) #background=numpy.array([0.00346008,0.00378121,0.00447618,0.00504666,0.00594937,0.00712394,0.00812853,0.00881078,0.0100297,0.0109190,0.0124697,0.0139781,0.0161290]) @@ -1442,7 +1465,7 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea #angle = float(key.split('_')[2].replace('Cos','')) half = key.split('_')[2].replace('Cos','') - angle = float(half.replace('.inc1.id1.sim','')) + angle = float(half.replace('.inc1.id1.sim','').replace(".gz", "")) if angle == angleSelection: # Get the number of reconstructed events @@ -1469,9 +1492,15 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea numberOfReconstructedEvents_pair = float(data[key][4][-2]) # Get the angular resolution - fwhm_tracked = data[key][2][7] - fwhm_untracked = data[key][3][7] - containment68 = data[key][4][6] + if data_nominal is not None: + fwhm_tracked = data_nominal[key][2][7] + fwhm_untracked = data_nominal[key][3][7] + containment68 = data_nominal[key][4][6] + else: + fwhm_tracked = data[key][2][7] + fwhm_untracked = data[key][3][7] + containment68 = data[key][4][6] + # Calculate the effective area effectiveArea_tracked = (numberOfReconstructedEvents_tracked/numberOfSimulatedEvents) * math.pi * SurroundingSphere**2 @@ -1646,7 +1675,7 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea else: mev2erg=1 unit='MeV' - ylim=[1e-7,5e-5] + ylim=[1e-8,1e-2] plot.scatter(Energy, Sensitivity_tracked*mev2erg, color='darkgreen') plot.plot(Energy, Sensitivity_tracked*mev2erg, color='darkgreen', alpha=0.5, label='Compton') @@ -1655,7 +1684,7 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea plot.scatter(Energy, Sensitivity_pair*mev2erg, color='darkred') plot.plot(Energy, Sensitivity_pair*mev2erg, color='darkred', alpha=0.5, label='Pair') plot.xlabel(r'Energy (MeV)') - plot.ylabel(r'$E^2 \times$ Intensity ('+unit+' cm$^{-2}$ s$^{-1}$ sr$^{-1}$)') + plot.ylabel(r'$E^2 \times$ Intensity ('+unit+' cm$^{-2}$ s$^{-1}$)') plot.ylim(ylim) if xlog: @@ -1667,11 +1696,38 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea plot.legend(numpoints=1, scatterpoints=1, fontsize=16, frameon=True, loc='upper left') if save: - plot.savefig('SourceSensitivity.pdf', bbox_inches='tight') - plot.savefig('SourceSensitivity.png', bbox_inches='tight') + plot.savefig('SourceSensitivity_{}.pdf'.format(angleSelection), bbox_inches='tight') + plot.savefig('SourceSensitivity_{}.png'.format(angleSelection), bbox_inches='tight') if doplot: plot.show() + + + # writing txt files + results_txt_TC = open( '%s_Sensi_%.1fs_Cos%s_TC.txt' % (txtOutfileLabel, exposure, angleSelection), 'w') + results_txt_UC = open( '%s_Sensi_%.1fs_Cos%s_UC.txt' % (txtOutfileLabel, exposure, angleSelection), 'w') + results_txt_P = open( '%s_Sensi_%.1fs_Cos%s_P.txt' % (txtOutfileLabel, exposure, angleSelection), 'w') + + unit = "erg" if uniterg else "MeV" + + results_txt_TC.write("# Energy[MeV] E2Sensitivity[%scm-2s-1]\n" % unit ) + for ii, en in enumerate(Energy): + results_txt_TC.write("%.3f\t%.5g\n"%(en, Sensitivity_tracked[ii])) + results_txt_TC.close() + print('Created %s_Sensi_%.1fs_Cos%s_TC.txt' % (txtOutfileLabel, exposure, angleSelection)) + + results_txt_UC.write("# Energy[MeV] E2Sensitivity[%scm-2s-1]\n" % unit ) + for ii, en in enumerate(Energy): + results_txt_UC.write("%.3f\t%.5g\n"%(en, Sensitivity_untracked[ii])) + results_txt_UC.close() + print('Created %s_Sensi_%.1fs_Cos%s_UC.txt' % (txtOutfileLabel, exposure, angleSelection)) + + results_txt_P.write("# Energy[MeV] E2Sensitivity[%scm-2s-1]\n" % unit ) + for ii, en in enumerate(Energy): + results_txt_P.write("%.3f\t%.5g\n"%(en, Sensitivity_pair[ii])) + results_txt_P.close() + print( 'Created %s_Sensi_%.1fs_Cos%s_P.txt' % (txtOutfileLabel, exposure, angleSelection)) + return Energy, Sensitivity_tracked, Sensitivity_untracked, Sensitivity_pair @@ -1751,7 +1807,8 @@ def plotAllSourceSensitivities(data, angleSelection=0.8, plotIdeal=False, xlog=T plot.gca().set_yscale('log') plot.gca().set_xlim([1e-2,1e6]) plot.gca().set_ylim(ylim) - plot.gca().set_xlabel('Energy (MeV)', fontsize=16, fontname="Helvetica") + #plot.gca().set_xlabel('Energy (MeV)', fontsize=16, fontname="Helvetica") + plot.gca().set_xlabel('Energy (MeV)', fontsize=16 ) plot.gca().set_ylabel(r'$3\sigma$ Continuum Sensitivity $\times\ E^2$ ($\gamma$ '+unit+' cm$^{-2}$ s$^{-1}$)', fontsize=16) plot.annotate('Fermi-LAT', xy=(1e3,(1e-6*mev2erg)),xycoords='data',fontsize=16,color='magenta') @@ -1777,7 +1834,8 @@ def plotAllSourceSensitivities(data, angleSelection=0.8, plotIdeal=False, xlog=T #SPI ind=numpy.arange(20,46) plot.plot(energy[ind],sens[ind]*mev2erg,color='green',lw=2) - plot.annotate('SPI', xy=(6e-2,(1e-4*mev2erg)),xycoords='data',fontsize=16,color='green') + plot.annotate('SPI', xy=(2e-2,(1e-5*mev2erg)),xycoords='data',fontsize=16,color='green') + #COMPTEL comptel_energy=numpy.array([0.73295844,0.8483429,1.617075,5.057877,16.895761,29.717747]) @@ -1844,17 +1902,17 @@ def plotAllSourceSensitivities(data, angleSelection=0.8, plotIdeal=False, xlog=T #Combine with activation simulation numbers - amego_activation_interp = numpy.interp(compair_eng[(compair_eng>0.5) & (compair_eng<=10)],amego_activation_energy, amego_activation_sens) + ##HFamego_activation_interp = numpy.interp(compair_eng[(compair_eng>0.5) & (compair_eng<=10)],amego_activation_energy, amego_activation_sens) #plot.plot(compair_eng[(compair_eng>0.3) & (compair_eng<=10)], amego_activation_interp,color='blue', lw=3, linestyle='--') - for e in range(len(compair_eng[(compair_eng>0.5) & (compair_eng<=10)])): - combined[e+2]= (amego_activation_interp[e] + combined[e+2])/2 + ##HFfor e in range(len(compair_eng[(compair_eng>0.5) & (compair_eng<=10)])): + ##HF combined[e+2]= (amego_activation_interp[e] + combined[e+2])/2 #plot.plot(compair_eng,tracked,'r--',color='grey',lw=2) #plot.plot(compair_eng,pair,'r--',color='grey',lw=2) - plot.plot(compair_eng[(compair_eng>0.1) & (compair_eng<=1e3)],combined[(compair_eng>0.1) & (compair_eng<=1e3)]*mev2erg,color='black',lw=4) + plot.plot(compair_eng[(compair_eng>0.05) & (compair_eng<=1e3)],combined[(compair_eng>0.05) & (compair_eng<=1e3)]*mev2erg,color='black',lw=4) - print("Final AMEGO Numbers, energy:", compair_eng[(compair_eng>0.1) & (compair_eng<=1e3)]) - print("Final AMEGO Numbers, energy:", combined[(compair_eng>0.1) & (compair_eng<=1e3)]*mev2erg ) + ##HFprint("Final AMEGO Numbers, energy:", compair_eng[(compair_eng>0.1) & (compair_eng<=1e3)]) + print("Final AMEGO Numbers, energy:", combined[(compair_eng>0.05) & (compair_eng<=1e3)]*mev2erg ) #if plotIdeal: #plot.plot(compair_eng,tracked_ideal,'r:',color='black',lw=3) @@ -1867,16 +1925,21 @@ def plotAllSourceSensitivities(data, angleSelection=0.8, plotIdeal=False, xlog=T #plot.annotate('Previous ComPair', xy=(7e2,3e-6),xycoords='data',fontsize=14,color='red') #plot.plot(compair_eng,pair_idealang,'r-.',color='black',lw=3) - plot.annotate('AMEGO', xy=(1,(3e-7*mev2erg)),xycoords='data',fontsize=22) + #plot.annotate('AMEGO', xy=(1,(3e-7*mev2erg)),xycoords='data',fontsize=22) + plot.annotate('AMEGO-X (AstroPix)', xy=(0.2,(8e-7*mev2erg)),xycoords='data',fontsize=20) if save: plot.savefig('full_sensitivity_Cos%s.pdf' % angleSelection) plot.savefig('full_sensitivity_Cos%s.png' % angleSelection) - fout=open('AMEGO_sensitivity_Cos%s.dat' %angleSelection,"w") - fout.write('Energy (MeV)'+"\t"+ 'E2xSensitivity ('+unit+'cm-2 s-1)') - for i in range(len(compair_eng[(compair_eng>0.1) & (compair_eng<=1e3)])): - fout.write(str(compair_eng[(compair_eng>0.1) & (compair_eng<=1e3)][i])+"\t"+\ - str(combined[(compair_eng>0.1) & (compair_eng<=1e3)]*mev2erg)) + fout=open('AMEGO-X_sensitivity_Cos%s.dat' %angleSelection,"w") + fout.write('#Energy (MeV)'+"\t"+ 'E2xSensitivity ('+unit+'cm-2 s-1)\n') + for i in range(len(compair_eng[(compair_eng>0.05) & (compair_eng<=1e3)])): + + fout.write("%.3f\t%.5g\n"%(compair_eng[(compair_eng>0.05) & (compair_eng<=1e3)][i], combined[i]*mev2erg)) + #fout.write(str(compair_eng[(compair_eng>0.05) & (compair_eng<=1e3)][i])+"\t"+\ + #str(combined[(compair_eng>0.1) & (compair_eng<=1e3)]*mev2erg)) + # str(combined[i]*mev2erg)+"\n") + fout.close() if doplot == True: plot.show() diff --git a/Prep_cosima_pointsource.py b/Prep_cosima_pointsource.py new file mode 100644 index 0000000..6897ed6 --- /dev/null +++ b/Prep_cosima_pointsource.py @@ -0,0 +1,73 @@ +# Adapted from https://github.com/ComPair/python/blob/master/Prep_cosima.py + +# This code will generate source files and runCosima.sh + + +import numpy as np +import os + +import argparse + + +#define your energies and angles + +Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7,4,4.2,4.5,4.7,5,5.2,5.5,5.7,6] #all energies +#Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7] #compton only + +#theta_angles =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles +theta_angles =[0, 36.9] #limited angle selection +phi_angles=[0.0] + +Log_E=np.array(Log_E) +theta_angles=np.array(theta_angles) +phi_angles=np.array(phi_angles) + +OneBeam= 'FarFieldPointSource' + + +#gives the cosTheta array +def ang2cos(allAng): + return np.round(np.cos(np.deg2rad(allAng)),1) + +def logE2ene(allEne): + return (10**allEne).astype(int) + +energies=logE2ene(Log_E) +cos_ang =ang2cos(theta_angles) + + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("-g", help="Geometry file") + parser.add_argument("-o", default=os.getcwd(), help="Output file dir (Default: current directory)") + args = parser.parse_args() + + output_dir= args.o + geometry_file = args.g + + os.makedirs(output_dir,exist_ok=True) + + with open("{0}/runCosima.sh".format(output_dir), mode='w') as f: + for myene in energies: + for cosTh,ang in zip(cos_ang,theta_angles): + for phi in phi_angles: + + # this is to print all the parameters combinations + #print (geofile, OneBeam, myene/1000., cosTh, OneBeam, ang, myene) + + #simulate 1e6 events below 10 MeV (compton range) and 1e5 events >= 10 MeV (pair range) + NTriggers = 100000 if myene > 10e3 else 1000000 + + #this is just a long string, with all the raws of the .source file, and the energies/angles values + string= "# An example run for Cosima \n# This was created with the python wrapper --> create_source_file.py <--\n\nVersion 1 \nGeometry %s // Update this to your path \nCheckForOverlaps 1000 0.01 \nPhysicsListEM LivermorePol \n\nStoreCalibrate true\nStoreSimulationInfo true\nStoreOnlyEventsWithEnergyLoss true // Only relevant if no trigger criteria is given! \nDiscretizeHits true \nStoreSimulationInfoIonization false \n\nRun FFPS \nFFPS.FileName %s/%s_%.3fMeV_Cos%.1f_Phi%.1f \nFFPS.NTriggers %d \n\n\nFFPS.Source One \nOne.ParticleType 1 \nOne.Beam %s %.1f %.1f \nOne.Spectrum Mono %i\nOne.Flux 10000.0 "%(geometry_file, output_dir, OneBeam, myene/1000., cosTh, phi, NTriggers, OneBeam, ang, phi, myene) + source_file='%s/%s_%.3fMeV_Cos%.1f_Phi%.1f.source'%(output_dir, OneBeam,myene/1000.,cosTh, phi) + sf=open(source_file,'w') + sf.write(string) + sf.close() + + f.write("cosima -z -v 0 -s 120 {}\n".format(source_file)) + + print("Run this script: {0}/runCosima.sh".format(output_dir) ) + + diff --git a/Prep_cosima_ring.py b/Prep_cosima_ring.py new file mode 100644 index 0000000..2841828 --- /dev/null +++ b/Prep_cosima_ring.py @@ -0,0 +1,74 @@ +# Adapted from https://github.com/ComPair/python/blob/master/Prep_cosima.py + +# This code will generate source files and runCosima.sh + + +import numpy as np +import os + +import argparse + +Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7,4,4.2,4.5,4.7,5,5.2,5.5,5.7,6] #all energies +#Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7] #compton only + +#angles =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles +angles =[0, 36.9] #limited angle selection + +width = 1.0 # degrees + +Log_E=np.array(Log_E) +angles=np.array(angles) + +OneBeam= 'FarFieldAreaSource' + + +#gives the cosTheta array +def ang2cos(allAng): + return np.round(np.cos(np.deg2rad(allAng)),1) + +def logE2ene(allEne): + return (10**allEne).astype(int) + +energies=logE2ene(Log_E) +cos_ang =ang2cos(angles) + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("-g", help="Geometry file") + parser.add_argument("-o", default=os.getcwd(), help="Output file dir (Default: current directory)") + args = parser.parse_args() + + output_dir= args.o + geometry_file = args.g + + os.makedirs(output_dir,exist_ok=True) + + + with open("{0}/runCosima.sh".format(output_dir), mode='w') as f: + for myene in energies: + for cosTh,ang in zip(cos_ang,angles): + + # this is to print all the parameters combinations + #print (geofile, OneBeam, myene/1000., cosTh, OneBeam, ang, myene) + + minTheta = max( ang - width, 0 ) + maxTheta = min( ang + width, 180 ) + minPhi = 0.0 + maxPhi = 360.0 + + #simulate 1e6 events below 10 MeV (compton range) and 1e5 events >= 10 MeV (pair range) + NTriggers = 100000 if myene > 10e3 else 1000000 + + + #this is just a long string, with all the raws of the .source file, and the energies/angles values + string= "# An example run for Cosima \n# This was created with the python wrapper --> create_source_file.py <--\n\nVersion 1 \nGeometry %s // Update this to your path \nCheckForOverlaps 1000 0.01 \nPhysicsListEM LivermorePol \n\nStoreCalibrate true\nStoreSimulationInfo true\nStoreOnlyEventsWithEnergyLoss true // Only relevant if no trigger criteria is given! \nDiscretizeHits true \nStoreSimulationInfoIonization false \n\nRun FFPS \nFFPS.FileName %s/%s_%.3fMeV_Cos%.1f \nFFPS.NTriggers %d \n\n\nFFPS.Source One \nOne.ParticleType 1 \nOne.Beam %s %.1f %.1f %.1f %.1f \nOne.Spectrum Mono %i\nOne.Flux 10000.0 "%(geometry_file, output_dir, OneBeam, myene/1000., cosTh, NTriggers, OneBeam, minTheta, maxTheta, minPhi, maxPhi, myene) + source_file='%s/%s_%.3fMeV_Cos%.1f.source'%(output_dir, OneBeam,myene/1000.,cosTh) + sf=open(source_file,'w') + sf.write(string) + sf.close() + + f.write("cosima -z -v 0 -s 120 {}\n".format(source_file)) + + print("Run this script: {0}/runCosima.sh".format(output_dir) ) + + diff --git a/Prep_revan_pointsource.py b/Prep_revan_pointsource.py new file mode 100644 index 0000000..734f925 --- /dev/null +++ b/Prep_revan_pointsource.py @@ -0,0 +1,93 @@ +# Adapted from https://github.com/ComPair/python/blob/master/Prep_revan.py +# This code will generate soft links of the source files +# into the output directory and generate runRevan.sh file. +# Please provide a geometry file & the full paths to the source and output directories. +# (The symlink is needed to be able to run revan on the same input files for multiple +# revan config files/geometry files without overwriting the results, since revan +# always puts the output files in the same directory as the input, and there does not +# seem to be a way to specify the file name. + +import numpy as np +import math +import argparse +import os + + +#to be adjusted before running +Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7,4,4.2,4.5,4.7,5,5.2,5.5,5.7,6] #all energies +#Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7] #compton only + +#theta_angles =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles +theta_angles =[0, 36.9] #limited angle selection +phi_angles=[0.0] + + +OneBeam= 'FarFieldPointSource' + +Log_E=np.array(Log_E) +angles=np.array(angles) +phi_angles=np.array(phi_angles) + + +#gives the cosTheta array +def ang2cos(allAng): + return np.round(np.cos(np.deg2rad(allAng)),1) + +#in keV [316,501,1000,1585, ... ] +def logE2ene(allEne): + return (10**allEne).astype(int) + + +energies=logE2ene(Log_E) +cos_ang =ang2cos(angles) + +def create_config_file(source_dir, source_file, output_dir, geometry_file): + with open('{}/{}.revan.cfg'.format(output_dir, source_file), mode='w') as cfg: + with open('./{}'.format(base_file)) as base: + for line in base.readlines(): + if line.find('') != -1: + cfg.write('{}\n'.format(geometry_file)) + elif line.find('') != -1: + cfg.write('{}/{}.inc1.id1.sim.gz\n'.format(source_dir, source_file)) + else: + cfg.write(line) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("-g", help="Geometry file") + parser.add_argument("-f", default=os.getcwd(), help="Source file dir (Default: current directory)") + parser.add_argument("-o", default=os.getcwd(), help="Output file dir (Default: current directory); will create symbolic links to input files if output directory is different from input") + parser.add_argument("-b", default='revan_AMEGO_X.cfg', help="Base revan config file, default: revan_AMEGO_X.cfg") + args = parser.parse_args() + + source_dir= args.f + output_dir= args.o + geometry_file = args.g + base_file= args.b + + if not os.path.exists( output_dir): + os.makedirs( output_dir ) + + args = parser.parse_args() + with open("{}/runRevan.sh".format(output_dir), mode='w') as f: + for myene in energies: + for cosTh,ang in zip(cos_ang,theta_angles): + for phi in phi_angles: + source_file='%s_%.3fMeV_Cos%.1f_Phi%.1f'%(OneBeam,myene/1000.,cosTh, phi) + + create_config_file(source_dir, source_file, output_dir, geometry_file) + + simfile_in="{0}/{1}.inc1.id1.sim.gz".format(source_dir, source_file) + simfile_out="{0}/{1}.inc1.id1.sim.gz".format(output_dir, source_file) + + if not os.path.exists( simfile_out ): + if not os.path.exists( simfile_in ): + print (simfile_in) + assert os.path.exists( simfile_in ) + os.symlink(simfile_in, simfile_out) + + f.write("revan -a --oi -n -f {3} -g {2} -c {0}/{1}.revan.cfg >& {0}/log_{1}.txt \n".format(output_dir, source_file, geometry_file, simfile_out) ) + + print("Run this script: {}/runRevan.sh".format(output_dir)) + diff --git a/Prep_revan_ring.py b/Prep_revan_ring.py new file mode 100644 index 0000000..9e673c3 --- /dev/null +++ b/Prep_revan_ring.py @@ -0,0 +1,88 @@ +# Adapted from https://github.com/ComPair/python/blob/master/Prep_revan.py +# This code will generate soft links of the source files +# into the output directory and generate runRevan.sh file. +# Please provide a geometry file & the full paths to the source and output directories. +# (The symlink is needed to be able to run revan on the same input files for multiple +# revan config files/geometry files without overwriting the results, since revan +# always puts the output files in the same directory as the input, and there does not +# seem to be a way to specify the file name. + +import numpy as np +import math +import argparse +import os + +Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7,4,4.2,4.5,4.7,5,5.2,5.5,5.7,6] #all energies +#Log_E=[1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3,3.2,3.5,3.7] #compton only + +#angles =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles +angles =[0, 36.9] #limited angle selection + + +OneBeam= 'FarFieldAreaSource' + +Log_E=np.array(Log_E) +angles=np.array(angles) + + +#gives the cosTheta array +def ang2cos(allAng): + return np.round(np.cos(np.deg2rad(allAng)),1) + +#in keV [316,501,1000,1585, ... ] +def logE2ene(allEne): + return (10**allEne).astype(int) + + +energies=logE2ene(Log_E) +cos_ang =ang2cos(angles) + +def create_config_file(source_dir, source_file, output_dir, geometry_file): + with open('{}/{}.revan.cfg'.format(output_dir, source_file), mode='w') as cfg: + with open('./{}'.format(base_file)) as base: + for line in base.readlines(): + if line.find('') != -1: + cfg.write('{}\n'.format(geometry_file)) + elif line.find('') != -1: + cfg.write('{}/{}.inc1.id1.sim.gz\n'.format(source_dir, source_file)) + else: + cfg.write(line) + + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument("-g", help="Geometry file") + parser.add_argument("-f", default=os.getcwd(), help="Source file dir (Default: current directory)") + parser.add_argument("-o", default=os.getcwd(), help="Output file dir (Default: current directory); will create symbolic links to input files if output directory is different from input") + parser.add_argument("-b", default='revan_AMEGO_X.cfg', help="Base revan config file, default: revan_AMEGO_X.cfg") + args = parser.parse_args() + + source_dir= args.f + output_dir= args.o + geometry_file = args.g + base_file= args.b + + if not os.path.exists( output_dir): + os.makedirs( output_dir ) + + args = parser.parse_args() + with open("{}/runRevan.sh".format(output_dir), mode='w') as f: + for myene in energies: + for cosTh,ang in zip(cos_ang,angles): + source_file='%s_%.3fMeV_Cos%.1f'%(OneBeam,myene/1000.,cosTh) + + create_config_file(source_dir, source_file, output_dir, geometry_file) + + simfile_in="{0}/{1}.inc1.id1.sim.gz".format(source_dir, source_file) + simfile_out="{0}/{1}.inc1.id1.sim.gz".format(output_dir, source_file) + + if not os.path.exists( simfile_out ): + if not os.path.exists( simfile_in ): + print (simfile_in) + assert os.path.exists( simfile_in ) + os.symlink(simfile_in, simfile_out) + + f.write("revan -a --oi -n -f {3} -g {2} -c {0}/{1}.revan.cfg >& {0}/log_{1}.txt \n".format(output_dir, source_file, geometry_file, simfile_out) ) + + print("Run this script: {}/runRevan.sh".format(output_dir)) + diff --git a/README.md b/README.md index 66273d0..ce93e00 100644 --- a/README.md +++ b/README.md @@ -2,4 +2,43 @@ python scripts to manipulate the output of MEGAlib For a description of the tools ad how to make performance plots, see https://confluence.slac.stanford.edu/display/COMPair/AMEGO+Python+Tools +New: +Use `Prep_cosima_ring.py` and `Prep_revan_ring.py` to simulate FarFieldAreaSources, with the directions of the incoming photons randomly distributed on a ring centered around the boresight, with the ring radius given by the theta values and the ring width set to one. Edit the scripts to adjust theta values and energies if necessary. + +usage: + +``` +Prep_cosima_ring.py [-h] [-g G] [-o O] + +optional arguments: + -h, --help show this help message and exit + -g G Geometry file + -o O Output file dir (Default: current directory) +``` + +``` +Prep_revan_ring.py [-h] [-g G] [-f F] [-o O] [-b B] + +optional arguments: + -h, --help show this help message and exit + -g G Geometry file + -f F Source file dir (Default: current directory) + -o O Output file dir (Default: current directory); will create + symbolic links to input files if output directory is different + from input + -b B Base revan config file, default: revan_AMEGO_X.cfg +``` + + +Use `doAnalysis.py` for all-in-one analysis (using EventAnalysis and FigureOfMeritPlotter). + +``` +usage: doAnalysis.py [-h] [-d DIR] [-t TAG] + +optional arguments: + -h, --help show this help message and exit + -d DIR, --dir DIR Directory with .sim and .tra files. (Default: current + directory) + -t TAG, --tag TAG Tag for output txt files. (Default: xxx) +``` diff --git a/doAnalysis.py b/doAnalysis.py new file mode 100644 index 0000000..0226b9d --- /dev/null +++ b/doAnalysis.py @@ -0,0 +1,29 @@ +import EventAnalysis, FigureOfMeritPlotter, os +import argparse + +parser = argparse.ArgumentParser() +parser.add_argument("-d", "--dir", default=os.getcwd(), help="Directory with .sim and .tra files. (Default: current directory)") +parser.add_argument("-t", "--tag", default="xxx", help="Tag for output txt files. (Default: xxx)") +args = parser.parse_args() + + +#this is needed so that the plots will be created in the directory with the .sim and .tra files, not the directory with the scripts. +os.chdir(args.dir ) + +if not os.path.exists( "./TriggerEfficiency.txt" ): + + EventAnalysis.getTriggerEfficiency( directory = "./" ) + +EventAnalysis.performCompleteAnalysis( directory="./", showPlots = False ) + +data = FigureOfMeritPlotter.parseEventAnalysisLogs("./", triggerEfficiencyFilename="./TriggerEfficiency.txt") + +sel= [1.0, 0.8] + +FigureOfMeritPlotter.plotAngularResolution(data, angleSelections=sel, save=True, txtOutfileLabel=args.tag) +FigureOfMeritPlotter.plotEnergyResolution(data, angleSelections=sel, save=True, txtOutfileLabel=args.tag) +FigureOfMeritPlotter.plotEffectiveArea(data, angleSelections=sel, ylog=True, save=True, txtOutfileLabel=args.tag) + +for angle in sel: + FigureOfMeritPlotter.plotSourceSensitivity(data, angleSelection=angle, doplot=True, save=True, uniterg=False, exposure = 0.2 * 3 * 365.25 * 24 * 3600, txtOutfileLabel=args.tag) +