From 0351faad53efe80f8cd934e4311fb0632f9c0443 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 26 Apr 2021 14:42:46 -0400
Subject: [PATCH 01/30] Updated scripts for cosima and revan

---
 Prep_cosima_pointsource.py | 66 ++++++++++++++++++++++++++++
 Prep_cosima_ring.py        | 67 ++++++++++++++++++++++++++++
 Prep_revan_pointsource.py  | 90 ++++++++++++++++++++++++++++++++++++++
 Prep_revan_ring.py         | 85 +++++++++++++++++++++++++++++++++++
 4 files changed, 308 insertions(+)
 create mode 100644 Prep_cosima_pointsource.py
 create mode 100644 Prep_cosima_ring.py
 create mode 100644 Prep_revan_pointsource.py
 create mode 100644 Prep_revan_ring.py

diff --git a/Prep_cosima_pointsource.py b/Prep_cosima_pointsource.py
new file mode 100644
index 0000000..0d7f87d
--- /dev/null
+++ b/Prep_cosima_pointsource.py
@@ -0,0 +1,66 @@
+# 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]
+theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+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)
+         
+                #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             100000 \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, 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..897f9a2
--- /dev/null
+++ b/Prep_cosima_ring.py
@@ -0,0 +1,67 @@
+# 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]
+angles = [0, 25.8, 36.9, 45.6, 53.1, 60] #cos
+
+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
+            
+                #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             100000 \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, 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..06baa90
--- /dev/null
+++ b/Prep_revan_pointsource.py
@@ -0,0 +1,90 @@
+# 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]
+theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+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('<GeometryFileName>') != -1:
+                    cfg.write('<GeometryFileName>{}</GeometryFileName>\n'.format(geometry_file))
+                elif line.find('<DataFileName>') != -1:
+                    cfg.write('<DataFileName>{}/{}.inc1.id1.sim.gz</DataFileName>\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..a04c29d
--- /dev/null
+++ b/Prep_revan_ring.py
@@ -0,0 +1,85 @@
+# 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]
+angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+
+
+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('<GeometryFileName>') != -1:
+                    cfg.write('<GeometryFileName>{}</GeometryFileName>\n'.format(geometry_file))
+                elif line.find('<DataFileName>') != -1:
+                    cfg.write('<DataFileName>{}/{}.inc1.id1.sim.gz</DataFileName>\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))
+                

From 05525bba3322bb677b5f8e9458659515336d30d0 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 26 Apr 2021 14:49:26 -0400
Subject: [PATCH 02/30] Support for gzipped files

---
 EventAnalysis.py | 20 +++++++++++++++-----
 1 file changed, 15 insertions(+), 5 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index c25dfd3..12b4dca 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
@@ -370,7 +371,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) )
@@ -2022,7 +2029,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') + glob.glob( './*.tra.gz')
 
 
     # Check if the user supplied a single file vs a list of files
@@ -2216,7 +2223,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 +2248,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 +2329,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:    

From 249b5eaa768191868e8c1b6c34cc4a9713ff034d Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 26 Apr 2021 15:01:59 -0400
Subject: [PATCH 03/30] use OI lines for true source direction

---
 EventAnalysis.py | 44 ++++++++++++++++++++++++++++++++++++++------
 1 file changed, 38 insertions(+), 6 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 12b4dca..94c95cd 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -340,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
@@ -397,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() 
 
@@ -436,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
@@ -606,6 +623,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
@@ -717,6 +738,8 @@ 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
 
     # Print some event statistics
     print("\n\nStatistics of Event Selection")
@@ -1029,14 +1052,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 )
 
-        # 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)
+        if direction_source is None:
 
-        # Set the origin position of the original gamma-ray
-        position_source = [position_conversion[0]-dx, position_conversion[1], 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]
+
+            # 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]

From 9acc50a499e1ebdcf4c6122790826a3e3596bbbf Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 17 May 2021 17:33:50 -0400
Subject: [PATCH 04/30] Default: Only 2 angles

---
 Prep_cosima_pointsource.py | 3 ++-
 Prep_cosima_ring.py        | 3 ++-
 Prep_revan_pointsource.py  | 3 ++-
 Prep_revan_ring.py         | 3 ++-
 4 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/Prep_cosima_pointsource.py b/Prep_cosima_pointsource.py
index 0d7f87d..4aeea7a 100644
--- a/Prep_cosima_pointsource.py
+++ b/Prep_cosima_pointsource.py
@@ -11,7 +11,8 @@
 
 #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]
-theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+#theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+theta_angles  =[0, 36.9]
 phi_angles=[0.0]
 
 Log_E=np.array(Log_E)
diff --git a/Prep_cosima_ring.py b/Prep_cosima_ring.py
index 897f9a2..a417941 100644
--- a/Prep_cosima_ring.py
+++ b/Prep_cosima_ring.py
@@ -9,7 +9,8 @@
 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]
-angles = [0, 25.8, 36.9, 45.6, 53.1, 60] #cos
+#angles = [0, 25.8, 36.9, 45.6, 53.1, 60] #cos
+angles = [0, 36.9] #cos
 
 width = 1.0 # degrees
 
diff --git a/Prep_revan_pointsource.py b/Prep_revan_pointsource.py
index 06baa90..4f2e295 100644
--- a/Prep_revan_pointsource.py
+++ b/Prep_revan_pointsource.py
@@ -15,7 +15,8 @@
 
 #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]
-theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+#theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+theta_angles  =[0, 36.9]
 phi_angles=[0.0]
 
 
diff --git a/Prep_revan_ring.py b/Prep_revan_ring.py
index a04c29d..07eab92 100644
--- a/Prep_revan_ring.py
+++ b/Prep_revan_ring.py
@@ -13,7 +13,8 @@
 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]
-angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+#angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
+angles  =[0, 36.9]
 
 
 OneBeam= 'FarFieldAreaSource'

From 9b862b8589d8a0e2b066d8eaaf3b66567d591fd0 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 17 May 2021 17:38:34 -0400
Subject: [PATCH 05/30] Documentation

---
 README.md | 26 ++++++++++++++++++++++++++
 1 file changed, 26 insertions(+)

diff --git a/README.md b/README.md
index 66273d0..f040b1d 100644
--- a/README.md
+++ b/README.md
@@ -2,4 +2,30 @@ 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
+```

From a98ed03944b5502f4ec1f778ca5eac3866af785c Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 17 May 2021 18:54:54 -0400
Subject: [PATCH 06/30] Small bug fixes

---
 EventAnalysis.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 94c95cd..06a8bd3 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1085,9 +1085,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)
@@ -2061,7 +2058,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( './*.tra') + glob.glob( './*.tra.gz')
+        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
@@ -2110,6 +2107,8 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
             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

From 0edf155e073b40c611b45499aad93921d9a60a7d Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 17 May 2021 18:55:16 -0400
Subject: [PATCH 07/30] Increase precision in output files; skip non-existent
 log files

---
 FigureOfMeritPlotter.py | 31 ++++++++++++++++---------------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/FigureOfMeritPlotter.py b/FigureOfMeritPlotter.py
index 660e421..8620595 100644
--- a/FigureOfMeritPlotter.py
+++ b/FigureOfMeritPlotter.py
@@ -258,16 +258,10 @@ 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')
 
@@ -279,6 +273,13 @@ 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
@@ -432,19 +433,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))
 
@@ -728,19 +729,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))
         
@@ -1056,19 +1057,19 @@ 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))
         

From ea391e867cd5c89950a7045048f104dc2b3de590 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 21 May 2021 13:14:30 -0400
Subject: [PATCH 08/30] Change default range for ARM fit

---
 EventAnalysis.py | 28 +++++++++++++++-------------
 1 file changed, 15 insertions(+), 13 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 06a8bd3..bfeaa72 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -847,12 +847,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')
@@ -860,11 +860,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])
@@ -2115,8 +2115,10 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
 
         # 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
 
             print("--------- All Compton Events ---------")
             # Calculate the energy resolution for Compton events
@@ -2197,7 +2199,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
             FWHM_pairComptonEvents = numpy.nan
 
         # Open the results filename for writing
-        output_filename = filename.replace('.tra','.log')
+        output_filename = filename.replace('.tra','.log').replace(".gz", "")
         output = open(output_filename, 'w')
 
         # Write the results to disk

From 7695e97af4b4d210ce37fd0258d2b78b7e620e0b Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 21 May 2021 13:15:28 -0400
Subject: [PATCH 09/30] Default exposure for Sensitivity curve: 3 years with
 20% duty cycle; txt output for sensitivity curve

---
 FigureOfMeritPlotter.py | 53 +++++++++++++++++++++++++++++++----------
 1 file changed, 40 insertions(+), 13 deletions(-)

diff --git a/FigureOfMeritPlotter.py b/FigureOfMeritPlotter.py
index 8620595..5a3cdbe 100644
--- a/FigureOfMeritPlotter.py
+++ b/FigureOfMeritPlotter.py
@@ -263,7 +263,7 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal
         numberOfSimulatedEvents = lineContents[2]
 
         # Generate the log filename
-        analysisLog = directory + '/' + simulationName.replace('.sim', '.log')
+        analysisLog = directory + '/' + simulationName.replace('.sim', '.log').replace(".gz", "")
 
         try:
             if silent==False:
@@ -394,7 +394,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)
@@ -563,7 +563,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:
@@ -689,7 +689,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)
@@ -862,7 +862,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)
@@ -988,7 +988,7 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False
             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:
@@ -1171,7 +1171,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
 
@@ -1278,7 +1278,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:
@@ -1400,7 +1400,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:
@@ -1422,9 +1422,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'):
 
     #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])
@@ -1443,7 +1443,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
@@ -1656,7 +1656,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:
@@ -1673,6 +1673,33 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3.536*10**7, idea
 
     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 
 

From f0599722b0a1bbc1f47cbce47c173828d1d44e8b Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 21 May 2021 13:16:20 -0400
Subject: [PATCH 10/30] all-in-one analysis script

---
 doAnalysis.py | 29 +++++++++++++++++++++++++++++
 1 file changed, 29 insertions(+)
 create mode 100644 doAnalysis.py

diff --git a/doAnalysis.py b/doAnalysis.py
new file mode 100644
index 0000000..4827e19
--- /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, save=True, uniterg=False, exposure = 0.2 * 3 * 365.25 * 24 * 3600, txtOutfileLabel=args.tag)
+

From e58617a8af8eac7aafa95e21f22fb947217f2e88 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 21 May 2021 13:18:28 -0400
Subject: [PATCH 11/30] Update documentation

---
 README.md | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/README.md b/README.md
index f040b1d..ce93e00 100644
--- a/README.md
+++ b/README.md
@@ -29,3 +29,16 @@ optional arguments:
               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)
+```

From 5fe9ec143d60b1aa5dcdd2229bc82c4e7ef429eb Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 21 May 2021 13:49:02 -0400
Subject: [PATCH 12/30] Plot sensitivity curve

---
 FigureOfMeritPlotter.py | 6 +++---
 doAnalysis.py           | 2 +-
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/FigureOfMeritPlotter.py b/FigureOfMeritPlotter.py
index 5a3cdbe..a505b14 100644
--- a/FigureOfMeritPlotter.py
+++ b/FigureOfMeritPlotter.py
@@ -1647,7 +1647,7 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3*365.25*24*3600*
         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')
@@ -1668,8 +1668,8 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3*365.25*24*3600*
         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()
diff --git a/doAnalysis.py b/doAnalysis.py
index 4827e19..0226b9d 100644
--- a/doAnalysis.py
+++ b/doAnalysis.py
@@ -25,5 +25,5 @@
 FigureOfMeritPlotter.plotEffectiveArea(data, angleSelections=sel, ylog=True, save=True, txtOutfileLabel=args.tag)
 
 for angle in sel:
-    FigureOfMeritPlotter.plotSourceSensitivity(data, angleSelection=angle, save=True, uniterg=False, exposure = 0.2 * 3 * 365.25 * 24 * 3600, txtOutfileLabel=args.tag)
+    FigureOfMeritPlotter.plotSourceSensitivity(data, angleSelection=angle, doplot=True, save=True, uniterg=False, exposure = 0.2 * 3 * 365.25 * 24 * 3600, txtOutfileLabel=args.tag)
 

From b52d2d1a67ca946a6a8de3d52c44f10b384eaaf8 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 2 Jun 2021 12:41:42 -0400
Subject: [PATCH 13/30] Fix bug (source direction for Compton events)

---
 EventAnalysis.py | 14 ++++++++++----
 1 file changed, 10 insertions(+), 4 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index bfeaa72..1d91861 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -571,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)
@@ -825,14 +831,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)

From bbdb2f3bf75fdba3bb9ced0572e7dd8bd4a1f8fc Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Tue, 15 Jun 2021 17:21:36 -0400
Subject: [PATCH 14/30] Fit range for low energies

---
 EventAnalysis.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 1d91861..7727d40 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -2125,6 +2125,8 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
             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

From 2f60a6df100c59f038cc48461102dbcd1cd160c5 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Tue, 15 Jun 2021 17:22:11 -0400
Subject: [PATCH 15/30] More statistics at lower energies

---
 Prep_cosima_pointsource.py | 14 ++++++++++----
 Prep_cosima_ring.py        | 16 +++++++++++-----
 Prep_revan_pointsource.py  |  8 +++++---
 Prep_revan_ring.py         |  8 +++++---
 4 files changed, 31 insertions(+), 15 deletions(-)

diff --git a/Prep_cosima_pointsource.py b/Prep_cosima_pointsource.py
index 4aeea7a..6897ed6 100644
--- a/Prep_cosima_pointsource.py
+++ b/Prep_cosima_pointsource.py
@@ -10,9 +10,12 @@
 
 
 #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]
-#theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
-theta_angles  =[0, 36.9]
+
+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)
@@ -52,9 +55,12 @@ def logE2ene(allEne):
          
                 # 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             100000 \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, OneBeam, ang, phi, myene)
+                    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)
diff --git a/Prep_cosima_ring.py b/Prep_cosima_ring.py
index a417941..30ac21f 100644
--- a/Prep_cosima_ring.py
+++ b/Prep_cosima_ring.py
@@ -8,9 +8,11 @@
 
 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]
-#angles = [0, 25.8, 36.9, 45.6, 53.1, 60] #cos
-angles = [0, 36.9] #cos
+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
 
 width = 1.0 # degrees
 
@@ -53,9 +55,13 @@ def logE2ene(allEne):
                 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             100000 \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, OneBeam, minTheta, maxTheta, minPhi, maxPhi, myene)
+                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)
diff --git a/Prep_revan_pointsource.py b/Prep_revan_pointsource.py
index 4f2e295..734f925 100644
--- a/Prep_revan_pointsource.py
+++ b/Prep_revan_pointsource.py
@@ -14,9 +14,11 @@
 
 
 #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]
-#theta_angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
-theta_angles  =[0, 36.9]
+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]
 
 
diff --git a/Prep_revan_ring.py b/Prep_revan_ring.py
index 07eab92..3b4cf8b 100644
--- a/Prep_revan_ring.py
+++ b/Prep_revan_ring.py
@@ -12,9 +12,11 @@
 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]
-#angles  =[0, 25.8, 36.9, 45.6, 53.1, 60]
-angles  =[0, 36.9]
+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
 
 
 OneBeam= 'FarFieldAreaSource'

From f7078eebb6d5f7fa78bbd0e9024f99444543f896 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 16 Jun 2021 12:19:00 -0400
Subject: [PATCH 16/30] Fix variable name

---
 Prep_cosima_ring.py | 4 ++--
 Prep_revan_ring.py  | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/Prep_cosima_ring.py b/Prep_cosima_ring.py
index 30ac21f..2841828 100644
--- a/Prep_cosima_ring.py
+++ b/Prep_cosima_ring.py
@@ -11,8 +11,8 @@
 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
+#angles  =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles
+angles  =[0, 36.9] #limited angle selection
 
 width = 1.0 # degrees
 
diff --git a/Prep_revan_ring.py b/Prep_revan_ring.py
index 3b4cf8b..9e673c3 100644
--- a/Prep_revan_ring.py
+++ b/Prep_revan_ring.py
@@ -15,8 +15,8 @@
 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
+#angles  =[0, 25.8, 36.9, 45.6, 53.1, 60] #all angles
+angles  =[0, 36.9] #limited angle selection
 
 
 OneBeam= 'FarFieldAreaSource'

From 5d6d82729ba6f9810aff5ce78111b284daa29e70 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 21 Jul 2021 19:12:19 -0400
Subject: [PATCH 17/30] new script for sensitivity calculation

---
 BurstSensitivities.py | 180 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 180 insertions(+)
 create mode 100644 BurstSensitivities.py

diff --git a/BurstSensitivities.py b/BurstSensitivities.py
new file mode 100644
index 0000000..e5c47e3
--- /dev/null
+++ b/BurstSensitivities.py
@@ -0,0 +1,180 @@
+#!/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
+
+__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
+
+
+
+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
+
+    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*(Emax - Emin), 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" )

From 424daea48b0c98ff9cd4822505ced2962b2f92dd Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Thu, 22 Jul 2021 09:13:26 -0400
Subject: [PATCH 18/30] Removed energy term

---
 BurstSensitivities.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/BurstSensitivities.py b/BurstSensitivities.py
index e5c47e3..9722e81 100644
--- a/BurstSensitivities.py
+++ b/BurstSensitivities.py
@@ -173,7 +173,7 @@ def get_sensitivity(time, significance, Aeff, BGrate):
 
         #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*(Emax - Emin), BGrate)
+            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, "

From 8c4dd8fc5ad0539c2f2be01dcab0ca93b87f1e48 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Fri, 20 Aug 2021 17:37:08 -0400
Subject: [PATCH 19/30] Include combination of all available channels

---
 BurstSensitivities.py | 43 ++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 42 insertions(+), 1 deletion(-)

diff --git a/BurstSensitivities.py b/BurstSensitivities.py
index 9722e81..a291587 100644
--- a/BurstSensitivities.py
+++ b/BurstSensitivities.py
@@ -54,7 +54,7 @@
 import matplotlib.pyplot as plt
 import ROOT
 from scipy.interpolate import interp1d
-from scipy.optimize import minimize
+from scipy.optimize import minimize, root_scalar
 
 __description__ = 'Integral flux sensitivities for bursts'
 
@@ -153,6 +153,15 @@ def get_sensitivity(time, significance, Aeff, BGrate):
     num = (significance**2/2+arg)/(Aeff*time)
     return num
 
+def get_signficance(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__':
@@ -165,6 +174,9 @@ def get_sensitivity(time, significance, Aeff, BGrate):
     
     Emin = args.emin
     Emax = args.emax
+    
+    EAs_for_combo = {}
+    BGs_for_combo = {}
 
     for type in theEAs.keys():
 
@@ -178,3 +190,32 @@ def get_sensitivity(time, 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_signficance( 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_signficance( 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 )

From 504a0e51e7caa823f7de73f11d7651ca5e6bc719 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 25 Aug 2021 16:49:04 -0400
Subject: [PATCH 20/30] typo

---
 BurstSensitivities.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/BurstSensitivities.py b/BurstSensitivities.py
index a291587..8edf937 100644
--- a/BurstSensitivities.py
+++ b/BurstSensitivities.py
@@ -153,7 +153,7 @@ def get_sensitivity(time, significance, Aeff, BGrate):
     num = (significance**2/2+arg)/(Aeff*time)
     return num
 
-def get_signficance(flux, time, Aeff, BGrate):
+def get_significance(flux, time, Aeff, BGrate):
     """Calculate expected significance for a given flux (ph/cm2/s), exposure time,
         effective area, and background rate"""
         
@@ -198,7 +198,7 @@ def get_combined_significance(flux):
     
         sig = 0
         for type in EAs_for_combo.keys():
-            sig += get_signficance( flux, args.exposure, EAs_for_combo[type], BGs_for_combo[type])**2
+            sig += get_significance( flux, args.exposure, EAs_for_combo[type], BGs_for_combo[type])**2
             
         return np.sqrt(sig)
         
@@ -211,7 +211,7 @@ def function_to_root( flux ):
     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_signficance( min_flux, args.exposure, EAs_for_combo[type], BGs_for_combo[type])
+        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]

From a3c5b3d7f873ad8c3cf980af253ea545b1f597f8 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 25 Aug 2021 16:49:19 -0400
Subject: [PATCH 21/30] Single-site evetns

---
 EventAnalysis.py        | 10 ++++--
 FigureOfMeritPlotter.py | 79 +++++++++++++++++++++++++++++------------
 2 files changed, 64 insertions(+), 25 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 7727d40..dffbac8 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -746,6 +746,7 @@ def parse(filename, sourceTheta=1.0, testnum=-1):
     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")
@@ -757,13 +758,15 @@ def parse(filename, sourceTheta=1.0, testnum=-1):
         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("")
 
@@ -2233,6 +2236,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
         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()
diff --git a/FigureOfMeritPlotter.py b/FigureOfMeritPlotter.py
index a505b14..dbc477d 100644
--- a/FigureOfMeritPlotter.py
+++ b/FigureOfMeritPlotter.py
@@ -283,6 +283,8 @@ def parseEventAnalysisLogs(directory, triggerEfficiencyFilename=None, silent=Fal
 
         # Loop through the analysis log file
         numberOfPairEventsIdeal = None
+        numberOfPhotoElectricEffectEvents = 0
+        
         for analysisLogLine in analysisLogLines:
 
             if "Results for simulation:" in analysisLogLine:
@@ -331,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
         
@@ -340,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:
@@ -958,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]
@@ -978,11 +983,13 @@ 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',''))
@@ -995,6 +1002,8 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False
                 #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
@@ -1022,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
@@ -1035,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)]
@@ -1042,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, \
@@ -1066,12 +1079,18 @@ def plotEffectiveArea(data, angleSelections=[1,0.9,0.8,0.7,0.6,0.5], ideal=False
         	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("%.3f\t%.3f\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
@@ -1088,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')
 
@@ -1424,7 +1446,7 @@ def resultsToFits(data, outfile='output.fits'):
 
 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, txtOutfileLabel='xxx'):
+    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])
@@ -1470,9 +1492,15 @@ def plotSourceSensitivity(data, angleSelection=0.8, exposure = 3*365.25*24*3600*
                 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
@@ -1779,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')
 
@@ -1805,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])
@@ -1872,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)
@@ -1895,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()

From b8e151e88bacd0bcedfe11c5131115d59d0291b9 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 25 Aug 2021 16:59:33 -0400
Subject: [PATCH 22/30] Extra case for PE/PP tokens

---
 EventAnalysis.py | 40 +++++++++++++++++++++-------------------
 1 file changed, 21 insertions(+), 19 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index dffbac8..59eead8 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -641,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:

From 6cd9c466860d27d8506c81ea54ee5da8c0cf9553 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Wed, 25 Aug 2021 17:01:08 -0400
Subject: [PATCH 23/30] Count single-site events

---
 EventAnalysis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 59eead8..4959030 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -755,7 +755,7 @@ def parse(filename, sourceTheta=1.0, testnum=-1):
     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

From 0f8f2d01ec2d551ba49ce12114ae96d90fa034cd Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 14:48:29 -0400
Subject: [PATCH 24/30] energy quantiles

---
 EventAnalysis.py | 21 ++++++++++++++++++---
 1 file changed, 18 insertions(+), 3 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 4959030..d1d7f92 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1513,7 +1513,7 @@ 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)")
@@ -1562,6 +1562,14 @@ 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))
+    
+    low, med, up = numpy.quantile( energy_pairReconstructedPhoton[selection], ( .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)
+
+    plot.savefig(fileBase+"_energyResolution_pairs.png")
 
     # Show the plot
     if showPlots == True:
@@ -1569,7 +1577,7 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N
     else:
         plot.close()
 
-    return fitMax, FWHM, sigma
+    return fitMax, FWHM, sigma, low, med, up
 
 ##########################################################################################
 
@@ -2199,7 +2207,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...")
@@ -2207,6 +2215,10 @@ 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
@@ -2235,6 +2247,9 @@ 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'])

From 19b96ef4059430732307936cd5d86274dfbd5afd Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 14:56:28 -0400
Subject: [PATCH 25/30] typo

---
 EventAnalysis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index d1d7f92..ec9e30e 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -2249,7 +2249,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
         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 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'])

From 112d42efccb6079556ec251c1ea096a14c595daa Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 15:17:04 -0400
Subject: [PATCH 26/30] typo

---
 EventAnalysis.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index ec9e30e..4b1c405 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1467,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
     
@@ -2132,6 +2132,10 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
         # Calculate the source theta in degrees
         source_theta = numpy.arccos(angle)*180./numpy.pi
 
+        output_filename = filename.replace('.tra','.log').replace(".gz", "")
+        if os.path.exists(output_filename):
+            continue
+
         # Don't bother measuring the energy and angular resolutuon values for Compton events above the specified maximumComptonEnergy
         if energy <= maximumComptonEnergy:
             
@@ -2224,7 +2228,6 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
             FWHM_pairComptonEvents = numpy.nan
 
         # Open the results filename for writing
-        output_filename = filename.replace('.tra','.log').replace(".gz", "")
         output = open(output_filename, 'w')
 
         # Write the results to disk

From f90a8fd02c732a9e16eb9802d0f057489666ec72 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 15:20:53 -0400
Subject: [PATCH 27/30] typo

---
 EventAnalysis.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 4b1c405..c26ab76 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -2122,6 +2122,10 @@ 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):
+            continue
+            
         if parsing:
             print("Parsing: %s %s Cos %s %s" % (energy, energySearchUnit, angle, filename))
             # Parse the .tra file obtained from revan
@@ -2132,9 +2136,6 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, angles
         # Calculate the source theta in degrees
         source_theta = numpy.arccos(angle)*180./numpy.pi
 
-        output_filename = filename.replace('.tra','.log').replace(".gz", "")
-        if os.path.exists(output_filename):
-            continue
 
         # Don't bother measuring the energy and angular resolutuon values for Compton events above the specified maximumComptonEnergy
         if energy <= maximumComptonEnergy:

From d9846f2dd767367c5ca52d9d8dab20b52b439de9 Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 15:24:51 -0400
Subject: [PATCH 28/30] fix return values

---
 EventAnalysis.py | 27 ++++++++++++++++-----------
 1 file changed, 16 insertions(+), 11 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index c26ab76..f8b90f1 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1502,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)]
@@ -1515,16 +1515,6 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N
     # Approximate sigma as 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
@@ -1577,6 +1567,20 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N
     else:
         plot.close()
 
+    # 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 - % 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
 
 ##########################################################################################
@@ -2124,6 +2128,7 @@ def performCompleteAnalysis(filename=None, directory=None, energies=None, 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:

From b099bf7ff214d97ba6ed1ab11cc838d9ec77d3ef Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 15:26:07 -0400
Subject: [PATCH 29/30] typo

---
 EventAnalysis.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index f8b90f1..0fbe5c4 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1574,7 +1574,7 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N
     print("")
     print("Fitting in range: ", energyFitRange[0], energyFitRange[1])
     print("Median reco energy: %s keV" % med)
-    print("68%% Containment Interval: %s - % keV" % (low, up) )
+    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)

From 114f88f241de403d2d6856ee0a76900edc2fb35a Mon Sep 17 00:00:00 2001
From: Henrike Fleischhack <henrike.fleischhack@googlemail.com>
Date: Mon, 1 Nov 2021 17:04:40 -0400
Subject: [PATCH 30/30] Fix energy selection

---
 EventAnalysis.py | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/EventAnalysis.py b/EventAnalysis.py
index 0fbe5c4..8870fe8 100644
--- a/EventAnalysis.py
+++ b/EventAnalysis.py
@@ -1536,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
@@ -1553,11 +1550,16 @@ def getEnergyResolutionForPairEvents(events, numberOfBins=100, energyPlotRange=N
     ax1.yaxis.set_minor_locator(AutoMinorLocator(4))
     ax2.xaxis.set_minor_locator(AutoMinorLocator(4))
     
-    low, med, up = numpy.quantile( energy_pairReconstructedPhoton[selection], ( .16, 0.5, .84) )
+    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")