fix errors

CHAOS.py

@@ -10,39 +10,33 @@ parser = argparse.ArgumentParser(prog = 'CHAOS', description = 'This program tra
 
 parser.add_argument('-s', '--size', default = 10, type = int, help = 'size of the simulation(number of rays)')
 
+#parser.add_argument('--xy', type = str, help = 'x and y argument for rastering along z-axis')
+
 parser.add_argument('--origin', type = str, help = 'origin of the bunch of rays')
 
 parser.add_argument('--height', default = 10, type = int, help = 'height of the scintillator')
 
 parser.add_argument('-e', '--edge', default = 51.962, type = float, help = 'edge length of the scintillator')
 
-parser.add_argument('--stop', default = 10000, type = float, help = 'cut off path length')
+parser.add_argument('--stop', default = 5000, type = float, help = 'cut off path length')
 
 parser.add_argument('-i', '--refraction_index', default = 2.15, type = float, help = 'index of refraction of the BGO cristall')
 
 parser.add_argument('-a', '--all', default = False, type = bool, help = 'parameter for rastering all')
 
-parser.add_argument('-x', '--x_value', default = False, type = int, help = 'vertical slice for rasterizing (from -25 to 25)')
+parser.add_argument('-x', '--x_value', default = False, type = str, help = 'vertical slice for rasterizing (from -25 to 25)')
 
 parser.add_argument('-o', '--outfn', default = None, help = 'Output file name')
 
 args = parser.parse_args()
 
+
 def check_range_x(value):
     ivalue = int(value)
     if ivalue < -50 or ivalue > 51:  # Modify the range as needed
         raise argparse.ArgumentTypeError(f"{value} is not in the valid range [-50,51]") 
     return ivalue
 
-if args.origin:
-    try:
-        x, y, z = map(float, args.origin.split(','))
-        origin = vector(x,y,z)
-    except ValueError:
-        print("Error: Invalid input format. Please provide the origin in the format 'x,y,z'")
-else: 
-    origin = vector(0,0,0)
-
 size = args.size
 height = args.height
 edge = args.edge
@@ -53,6 +47,14 @@ x_value = args.x_value
 check_range_x(x_value)
 outfn = args.outfn
 
+if args.origin:
+    try:
+        x, y, z = map(float, args.origin.split(','))
+        origin = vector(x,y,z)
+    except ValueError:
+        print("Error: Invalid input format. Please provide the origin in the format 'x,y,z'")
+
+
 if all:
     pos = 0
     origin = np.empty((145992,3,1))
@@ -85,27 +87,29 @@ if all:
             pos = pos + o
 
 if x_value:
-    if x_value < 25 and x_value > -25:
+    x = int(x_value)
+    if x < 25 and x > -25:
         apothem = int(np.sin(np.radians(60))*edge)
-        origin = np.empty((apothem*2*(height*2+1),3,1))
-        origin[:,0] = x_value
-        k = np.empty((apothem*2,1))
-        k[:,0] = np.array(np.arange(-apothem, apothem, 1))
-        for i in range(height*2+1):        
-            origin[apothem*2*i:(i+1)*apothem*2,1] = k
-            origin[apothem*2*i:(i+1)*apothem*2,2] = height - i
+        origin = np.empty(((apothem*2-1)*(height*2),3,1))
+        origin[:,0] = x
+        k = np.empty((apothem*2-1,1))
+        k[:,0] = np.array(np.arange(-apothem+0.5, apothem-0.5, 1))
+        for i in range(height*2):      
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,1] = k
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,2] = height - 0.5 - i
     else:
-        if x_value >= 25:
-            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(x_value-edge/2)))
+        if x >= 25:
+            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(x-edge/2)))
         else:
-            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*x_value))
-        origin = np.empty((apothem*2*(height*2+1),3,1))
-        origin[:,0] = x_value
-        k = np.empty((apothem*2,1))
-        k[:,0] = np.array(np.arange(-apothem, apothem, 1))
-        for i in range(height*2+1):
-            origin[apothem*2*i:(i+1)*apothem*2,1] = k
-            origin[apothem*2*i:(i+1)*apothem*2,2] = height - i
+            apothem = int(np.sin(np.radians(60))*edge - abs(np.tan(np.radians(60))*(-edge/2-x)))
+        origin = np.empty(((apothem*2-1)*(height*2),3,1))
+        origin[:,0] = x
+        k = np.empty((apothem*2-1,1))
+        k[:,0] = np.array(np.arange(-apothem+0.5, apothem-0.5, 1))
+        for i in range(height*2):
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,1] = k
+            origin[(apothem*2*i)-i:((i+1)*apothem*2)-i-1,2] = height - 0.5 - i
+
 
 
 det=[Plane(vector(-10,-35,5),vector(-10,-35,-5),vector(-10,-30,0),[]),
@@ -122,25 +126,43 @@ if outfn is None or outfn=='-':
 else: 
     csvfile = open(outfn, "w" )
 
-    
-r = Rays(size, origin)
-#calculate the incident th and phi angles
-th, phi = Rays.angle(r)
-sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
-#create origin vector for output file
-n = np.empty((len(origin)*size,1,3))
-if len(origin)>3:
-    for i in range(len(origin)):
-        n[size*i:size*(i+1),0] = origin[i].T
-else: 
-    n[...,0] = origin[0]
-    n[...,1] = origin[1]
-    n[...,2] = origin[2]
-#edit hit_rec so it can be written in one row
-hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
-#create data list
-data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
-#write output in csv_file
-for row in data:
-    csvfile.write(','.join(map(str, row)) + '\n')
+if x_value:
+    x = int(x_value)
+    for i in range(int(len(origin)/(2*height))):
+        r = Rays(size, origin[2*height*i:2*height*(i+1)])
+        n = np.empty((2*height*size,1,3))
+        n[...,0] = r.origin[:,0]
+        n[...,1] = r.origin[:,1]
+        n[...,2] = r.origin[:,2]
+        #calculate the incident th and phi angles
+        th, phi = Rays.angle(r)
+        sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+        #edit hit_rec so it can be written in one row
+        hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+        #create data list
+        data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
+        #write output in csv_file
+        for row in data:
+            csvfile.write(','.join(map(str, row)) + '\n')
+
+else:
+    r = Rays(size, origin[2*height*i:2*height*(i+1)])
+    th, phi = Rays.angle(r)
+    sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
+    #create origin vector for output file
+    n = np.empty((len(origin)*size,1,3))
+    if len(origin)>3:
+        for i in range(len(origin)):
+            n[size*i:size*(i+1),0] = origin[i].T
+    else: 
+        n[...,0] = origin[0]
+        n[...,1] = origin[1]
+        n[...,2] = origin[2]
+    #edit hit_rec so it can be written in one row
+    hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
+    #create data list
+    data = list(zip(n, th, phi, sc, tr, path, hit_rec_processed))
+    #write output in csv_file
+    for row in data:
+        csvfile.write(','.join(map(str, row)) + '\n')
 

filter_not.py

@@ -34,6 +34,10 @@ def main(csv_file, stop, Ind_S, height, edge):
     new_dir[:,0]= (np.sin(np.radians(th)) * np.cos(np.radians(phi))).reshape(-1,1)
     new_dir[:,1] = (np.sin(np.radians(th)) * np.sin(np.radians(phi))).reshape(-1,1)
     new_dir[:,2] = (np.cos(np.radians(th))).reshape(-1,1)
+    new_o = np.empty((len(origin),3,1))
+    new_o[:,0:3] = origin[:,:,0:3]
+    new_o = np.transpose(new_o, (0,2,1))
+    print(new_o)
     r = Rays_custom(origin, new_dir)
     csvfile_new = sys.stdout
     det=[Plane(vector(-10,-35,5),vector(-10,-35,-5),vector(-10,-30,0),[]),
@@ -45,7 +49,7 @@ def main(csv_file, stop, Ind_S, height, edge):
     th, phi = Rays_custom.angle(r)
     sc, tr, path, hit_rec = run(p, r, stop, Ind_S)
     hit_rec_processed = [','.join(','.join(map(str, row)) for row in arr) for arr in hit_rec]
-    data = list(zip(np.transpose(origin, axes=(0, 2, 1)), th, phi, sc, tr, path, hit_rec_processed))
+    data = list(zip(new_o, th, phi, sc, tr, path, hit_rec_processed))
     for row in data:
         csvfile_new.write(','.join(map(str, row)) + '\n')